Simplified exponential approximation for thermodynamics of a hard-core repulsive Yukawa fluid

Exponential approximation based on the first order mean spherical approximation (FMSA) is applied to the study of the structure and thermodynamics of hard-core repulsive Yukawa fluids. The proposed theory utilizes an exponential enhancement of the analytical solution of the FMSA due to Tang and Lu [J. Chem. Phys., 1993, 99, 9828] for the radial distribution function. From comparison with computer simulation data we have shown that at low density and low temperature conditions, where original FMSA theory fails, the FMSA-based exponential theory predicts a significant improvement.


Introduction
Authors are pleased and honored to contribute to the festschrift dedicated to the sixtieth anniversary of Orest Pizio. Around thirty years ago one of us (AT) started his scientific career under Orest's supervision. At that time we were dealing primarily with the electroneutral mixture of charged hard spheres that are widely known as a restricted primitive model of electrolytes. Such a model incorporates both Coulomb repulsion and Coulomb attraction. Repulsive interaction in the form of Yukawa potential, that is the subject of the present study, is quite frequently referred to as the screened Coulomb repulsion, and is widely used in modelling various real substances in statistical physics of soft matter. In particular, one can apply the Yukawa repulsion to the modelling of electrostatic interactions in one-component plasmas [1], electrical double layer interactions in charge-stabilized colloidal dispersions [2][3][4] and micellar surfactant solutions [5], etc. A brief but rather comprehensive review on the application of a repulsive Yukawa model in statistical physics can be found in the paper by Hopkins et al [6]. Some of these applications concern a pure repulsive Yukawa model, but very often the Yukawa repulsion is only a part of the potential functions used in modelling a wide class of the so-called complex fluids. Among others we wish to mention the core-softened fluid models, both homogeneous and those confined to a single surface, slit pore or even porous media. The field of core-softened fluid models is intensively developing nowadays and represents one of those in which Orest Pizio and his colleagues have recently taken interest [7][8][9].
Most generally, the repulsive Yukawa model is defined by the following potential function where γ > 0 is the energy (coupling) parameter that characterizes the magnitude of repulsion while parameter z controls the decay of repulsive interaction. This model has been extensively studied by means of various methods of statistical physics, such as numerical solutions of Ornstein-Zernike integral equation [3,4,10], density functional theory [11,12] and computer simulations [1,2,10,[12][13][14][15][16]. In contrast to the attractive Yukawa potential, characterized by γ < 0, the repulsive Yukawa model in general does not require the introduction of an extra core repulsion at short separations. However, looking for the analytical theory of the repulsive Yukawa fluid, attention was turned to the Waisman's solution of the mean-spherical approximation (MSA) for Yukawa interaction [17]. For reasons of the MSA solution technique, the Yukawa potential is written in a form that slightly differs from that of equation (1), namely, where σ is the hard-core diameter, ε = γ exp(−zσ) corresponds to the potential energy at a hard-core contact distance, r = σ . Since the MSA solution is not sensitive to the sign of energy parameter ε (or γ ), it can be applicable to both repulsive and attractive Yukawa interactions. The MSA theory for Yukawa fluids being analytical in general, requires a solution of a set of rather complex algebraic equations. Although the original solution due to Waisman [17] was later on simplified by Blum and Hoye [18], Cummings and Smith [19], Ginoza [20] and finally by Henderson et al [21] and Scalise et al [22], it still invokes a quadratic algebraic equation from which a unique physical root corresponding to the coupling parameter (prototype of the Debye screening parameter) should be deduced. Recently, looking for further simplification of the MSA theory, Tang and Lu [23] have solved Ornstein-Zernike equation by means of perturbation theory and developed the so-called first-order mean spherical approximation (FMSA) theory. While the FMSA is a simpler theory, it is only slightly less accurate than the full MSA solution, and in contrast to the latter, is "solvable" at any thermodynamic conditions [24]. Its simplicity, however, makes it also attractive for an employment in different related approaches, such as density functional theory [25,26] or augmented perturbation theory [27], etc.
Similar to the MSA theory, the FMSA solution by Tang and Lu [23] is valid for Yukawa potential (2) in both ε < 0 and ε > 0 cases. However, in the present study our interest will lie in applying the FMSA theory to the case ε > 0 only, and we will refer to this model as the hard-core repulsive Yukawa (HCRY) fluid. It has been already shown [10,28] that FMSA theory is reasonably accurate for dense HCRY fluids and at high temperatures, i.e., when Yukawa repulsion is weak or dominated by a hard-core repulsion. However, the main difficulty of the MSA-like theories for repulsive potentials concerns the appearance of negative values of the radial distribution function at short separations in the case of the growing strength of repulsive interactions. Thus, although this has not been documented in literature, it is expected that FMSA theory will start to fail if temperature in the HCRY fluid becomes lower. In particular, in what follows we will show that such a deficiency of the FMSA theory becomes especially drastic at low temperature and low density conditions. It is obvious that in this region of density and temperature parameters the problems occur not only with the FMSA predictions for the radial distribution function but also for the thermodynamics of the HCRY fluid. A natural way of fixing such a flaw in any linear theory is to employ an exponential approximation for the radial distribution function [29][30][31][32]. Besides ensuring positive values of the radial distribution function at short separations, the exponential approximation also provides a correct asymptote of the radial distribution function at large distances. This significant improvement for the radial distribution function can be extended towards thermodynamics by means of the so-called energy route when the internal energy of the system is evaluated using its standard definition through the radial distribution function. The main purpose of the present study is to employ this strategy for the HCRY fluid based on the existing solution of the FMSA theory for Yukawa potential. In the following section 2 we will introduce and determine what we are referring to as the simplified exponential (SEXP) theory based on the FMSA solution due to Tang and Lu [23]. This will be done for both the radial distribution function and for the thermodynamics. The results and their discussions are the subject of section 3. Since the low temperature and low density states have not been studied for the HCRY fluid so far, to verify the performance of the SEXP/FMSA theory we were forced to carry out the corresponding computer simulations. Thus, the necessary details of the canonical and constant pressure ensembles Monte Carlo techniques are briefly outlined as well. Our concluding remarks are summarized in the last section 4. To keep the paper simple and easy to understand, some auxiliary equations are moved to the appendix.

Exponential theory based on the FMSA solution
The structure and thermodynamics of the HCRY fluid, defined according to equation (2), will be studied here within the exponential approximation based on the analytical solution of the first-order meanspherical approximation (FMSA). The complete scheme and basic expressions of the FMSA solution for the Yukawa fluid are presented in separate papers by Tang and Lu [23,28]. Here we recall only the equations necessary for the radial distribution function.

FMSA solution for the radial distribution function
Like in any linear theory, the radial distribution function within the FMSA theory is represented as a sum of two terms, where β = (k B T ) −1 is the inverse of temperature, while g 0 (r ) is the radial distribution function of the reference system that is a fluid of pure hard spheres, and g 1 (r ) is the distance dependent correction term. The FMSA theory due to Tang and Lu provides us with the Laplace transform G 1 (s) for the correction term [23,28] where function Q 0 (s) refers to the Laplace transform of the Baxter factorization function resulting from the analytical solution of the Percus-Yevick (PY) approximation for pure hard spheres [34], with η = πρσ 3 /6 being the packing fraction parameter, while L(s) and S(s) are polynomials. The equation for the Laplace transform, G 0 (s), of the radial distribution function of a fluid of pure hard spheres within the framework of the PY reads [33,34] G 0 (s) For the purpose of the present study it is important to have rather compact expressions in order to invert the two functions, G 1 (s) and G 0 (s), into a real space. In the case of pure hard spheres, it reads [35,36] g 0 (r ) A similar inversion equation for the FMSA correction term is given by [32,36] The functions C (n 1 , n 2 , n 3 , r /σ) and D (n 1 , n 2 , n 3 , zσ, r /σ) as well as the still undefined polynomials L(s) and S(s) are introduced in appendix.

SEXP/FMSA approximation for the radial distribution function
Equations (3)-(8) represent the FMSA theory for the radial distribution function of the hard-core Yukawa fluid in general, i.e., for both the repulsive and attractive types of interaction. When evaluated numerically, the results predicted by equations (3)-(8) are both qualitatively and quantitatively quite similar to those obtained within the complete MSA theory and are rather satisfactory in all practical applications in the case of the attractive Yukawa interaction, i.e., when ε < 0. However, the same FMSA theory as well as the original MSA theory for the radial distribution function are not always that efficient in the case of repulsive Yukawa interaction, i.e., for the HCRY fluid which is the subject of the present study. In particular, by lowering the temperature for the HCRY fluid, one is actually increasing the strength of repulsive interaction. For such conditions, the radial distribution function g (r ) , being evaluated by means of equations (3)- (8), becomes negative at short distances r between particles. This is illustrated in figures 1, 2 and 3, where the results obtained within the FMSA theory are marked by dashed lines. Obviously, these physically incorrect results for the radial distribution function will also affect the FMSA predictions for the thermodynamics of the HCRY fluid. A natural way of fixing such a drawback of the linear theory is to employ an exponential (EXP) approximation for the radial distribution function as it has been suggested within the MSA theory [31]. In the case of the FMSA theory, this approximation reads g SEXP (r ) = g 0 (r )exp −βεg 1 (r ) . (9) Since FMSA theory itself is a simplified version of the more general MSA theory, the exponential approximation based on the FMSA solution is usually referred to as the simplified exponential (SEXP) approximation [29]. By comparing equations (9) and (3) one can immediately conclude that g SEXP (r ) is always positive, which is due to the FMSA correction, βεg 1 (r ) , residing in the exponent. These results are shown by solid lines in figures 1, 2 and 3. Additionally, exponential approximation (9) results in the radial distribution function that is asymptotically correct in the low density limit as well. Obviously, all these improved features of the exponential approximation remain valid within the original MSA theory. The reason for considering the SEXP/FMSA theory herein is due to the simplicity of the FMSA solution at nearly the same level of accuracy as the original MSA solution. This circumstance has caused a growing popularity of the FMSA theory itself. To illustrate this claim we are pointing out that in contrast to the MSA theory, the correction function g 1 (r ) given by equation (8) does not depend on the temperature. Hence, there is very simple linear temperature dependence of the FMSA correction for the radial distribution function. This feature will be exploited herein below to write down the thermodynamics within the SEXP/FMSA theory.  23003-5

Thermodynamics within the SEXP/FMSA approximation
Thermodynamics within the SEXP/FMSA theory will be obtained herein through the so-called energy route [30]. Following this scheme, the derivation of thermodynamics starts with the evaluation of the internal energy E using its definition through the radial distribution function, This result is used to obtain the Helmholtz free energy in the form where A 0 is the Helmholtz free energy of the fluid of hard spheres. Since functions g 0 (r ) and g 1 (r ) are independent of the temperature (see equation (9) for g SEXP (r )), the inverse temperature integral in the second term of equation (11) can be evaluated analytically, where we introduced notation ϕ 0 (x) = 1 − exp(−x) /x. Then, it is more convenient to rewrite the resulting Helmholtz free energy expression in the form where we have introduced the conventional FMSA Helmholtz free energy contribution with coefficients a 0 , a 1 and a 2 for the HCRY fluid given by [28,32] and the correction due to the exponential approximation is In contrast to equation (14), the integral in equation (18) should be evaluated numerically by integrating along the radial distance, r , using the Gaussian quadratures. The pressure, P , can be computed using a standard thermodynamic relation, where µ is the chemical potential. Within the SEXP/FMSA theory, the chemical potential can be represented in the form, This expression for chemical potential of the HCRY fluid is composed of two parts. The first part refers to the chemical potential within the FMSA theory, while the second one corresponds to the correction due to an exponential approximation for the radial distribution function. These two terms can be obtained from the following equations, where ϕ 1 (x) = 1 − x − exp (−x) /x 2 , and five density derivatives were introduced, i.e., ∂a 0 /∂ρ, ∂a 1 /∂ρ, ∂a 2 /∂ρ, ∂g 0 (r )/∂ρ and ∂g 1 (r )/∂ρ . The first two of them could be straightforwardly evaluated using definitions (15)-(17) for the quantities a 0 , a 1 and a 2 , respectively. As for the next two, in accordance with equations (7) and (8), one will be required to evaluate the density derivatives of the functions C (n 1 , n 2 , n 3 , x) and D n 1 , n 2 , n 3 , y, x . These derivatives are also straightforward and the details can be found in the appendix.

Results and discussions
We considered three models for the HCRY fluids with a fixed strength of repulsive energy, ε = 1, but three different ranges of repulsive interaction, i.e., characterized by three different values of the decay parameter, namely, zσ = 5, 3 and 1.8. The reduced temperature was fixed at T * = k B T /ε = 0.125. While three values of the decay parameter z have been already used in related studies [10,28] of the HCRY fluids, such a low value of reduced temperature has not been explored so far in theoretical studies of this type of model fluids. At the same time, the reduced temperature T * = 0.125 and even lower temperatures are quite obvious in computer simulation studies of the phase diagrams of the HCRY fluid [4,14]. While considering such a low temperature, our main intention, on the one hand, is to clearly illustrate the problems experienced by linear theories in the description of purely repulsive fluids, and, on the other hand, to show an extent to which the application of exponential theories can improve the theoretical treatment of repulsive fluids.
To test the theoretical predictions, we performed canonical ensemble and constant pressure ensemble Monte-Carlo (MC) simulations. To this end, the modified codes from the book of Frenkel and Smit [37] were used and the radial distribution function g (r ) , the compressibility factor Z = βP /ρ and the internal energy βE were calculated. In particular, the radial distribution function was obtained from canonical ensemble MC simulations for 256 particles in a cubic simulation box with periodic boundary conditions. All these simulations were performed for 500 equilibration steps and 1500 production steps. Each step consisted of 2500 attempts to displace a particle. On the other hand, constant pressure ensemble MC simulations were used to sample the density dependence of the pressure and internal energy. This set of simulations was performed for 256 particles with 2500 equilibration steps and 5000 production steps. Similarly, each step consisted of 10000 attempts to displace a particle and 100 attempts to change the box volume.
First, for each of three considered HCRY fluids we compare the radial distribution functions evaluated within the FMSA and SEXP/FMSA theories against MC simulations data. These results are shown in figures from 1 to 3 for the HCRY fluids characterized by zσ = 5, 3 and 1.8, respectively. Each part of these figures consists of a set of two theoretical (lines) and one computer simulation (symbols) results for one of four densities, namely, ρσ 3 = 0.1, 0.2, 0.4 and 0.6. As it was already mentioned, the main problem of the radial distribution functions evaluated within the linear theories, such as MSA and FMSA, concerns an appearance of the non-physical negative values of g (r ) at short distances r , when the strength of repulsion is increasing. Indeed, as we can see from figures 1, 2 and 3, this is the case at the considered reduced temperature T * = 0.125, where radial distribution functions shown by dashed lines (i.e., those that correspond to the FMSA theory) lie in the region of negative values. This region of negative values of g (r ) is more notable and extends far beyond the contact distance r /σ = 1 for lower densities ρσ 3 = 0.1 and 0.2. We note that for these two densities, the distances r for which g (r ) assume negative values, increase with an increase of the range of repulsion, i.e., going from zσ = 5, to zσ = 3 and 1.8. However, for density ρσ 3 = 0.4 and larger densities, the tendency changes to the opposite. Namely, the region of negative values of g (r ) of the HCRY fluid characterized by zσ = 5 (see figure 1), starts to shrink going to fluids with zσ = 3 and 1.8. In the case of the highest density explored in this study, ρσ 3 = 0.6, the radial distribution function g (r ) , evaluated from the FMSA theory, becomes even fully positive for zσ = 3 and 1.8 (see figures 2 and 3), while it remains always negative for the case of zσ = 5. More information concerning the radial distribution functions can be found in table 1, where the results for contact values, g (r = σ), are collected. Looking for the MC data one can see, that for each of three HCRY fluids the contact value g (r = σ) increases when density increases. However, in the case of the HCRY fluid with zσ = 5 the radial distribution function at the contact always remains smaller than one, g (r = σ) < 1, while in the case of two other HCRY fluids with zσ = 3 and 1.8 this is not the case. From the results presented in table 1 it follows that SEXP/FMSA theory results in g (r = σ) values that are always slightly lower than similar MC data, while in the case of the FMSA theory this depends on the density: at low densities FMSA theory results in contact values g (r = σ) that are lower than MC data, but at high densities it reverses, i.e., the MC data for g (r = σ) are smaller than those from the FMSA theory.
Thus, there is a crossover density which depends on the range of repulsion (i.e., depends on the value of parameter z ). The existence of this crossover density explains why for some HCRY fluid there is a range of densities, where the FMSA theory results in the contact values g (r = σ) that are very close to the MC data.
However, it is evident that on average in all considered cases the SEXP/FMSA approximation improves an agreement between theory and computer simulations, especially, at short distances. At intermediate and large distances, both FMSA and SEXP/FMSA theories are quite similar, and both slightly underestimate the ordering that takes place in the HCRY fluids. In particular, the oscillations of the radial distribution functions predicted by the FMSA and SEXP/FMSA theories show to be slightly weaker than those resulting from computer simulation data.
Corresponding results for thermodynamic properties are presented in figures 4, 5 and 6 and collected in tables 2, 3 and 4. In these figures and tables we show a comparison of the density dependence of the internal energy and compressibility factor for all three HCRY fluids, obtained from the FMSA and

Conclusions
Summarizing, a simplified exponential SEXP/FMSA theory was employed to describe the thermodynamics of the hard core repulsive Yukawa (HCRY) fluid. The theory is built up on an analytical solution of the FMSA theory for the radial distribution function due to Tang and Lu [23]. The main reason that has pursued us into this project were nonphysical negative values for the radial distribution function of the HCRY fluid that are predicted by the original FMSA theory at low temperatures. Obviously, this artifact of the FMSA theory will cause problems for the FMSA thermodynamics of the HCRY fluid as well.
The results presented in this study are obtained for the low reduced temperature characterized by T * = 0.125. Such a low reduced temperature for the HCRY fluids has not been explored within the integral equation theory studies so far. However, there are computer simulation studies of the phase diagram of the HCRY fluid at such and even lower temperatures [4,14]. We have shown that the conventional FMSA theory at low density leads to the non-physical negative values of the radial distribution function at short distances and, consequently, to incorrect thermodynamics yielding the negative internal energy, but it seems to be reasonably accurate at higher densities. By contrast, the exponential approximation based on the FMSA solution significantly improves the theoretical treatment of the HCRY fluids. The SEXP/FMSA theory proves to be more successful for the HCRY fluids characterized by a more short-ranged repulsive interaction. In the case of long-ranged repulsive potentials, the exponential approximation performing reasonably well at low densities, tends to predict weaker oscillations of the radial distribution function at higher densities.