Microscopic structure and thermodynamics of a core-softened model fluid from the second-order integral equations theory

We have studied the structure and thermodynamic properties of isotropic three-dimensional core-softened fluid by using the second-order Ornstein-Zernike integral equations completed by the hypernetted chain and Percus-Yevick closures. The radial distribution functions are compared with those from singlet integral equations and with computer simulation data. The limits of the region of density anomaly resulting from different approximate theories are established. The obtained results show that the second-order hypernetted chain approximation can be used to describe both the structure and the density anomaly of this model fluid. Moreover, we present the results of calculations of the bridge functions.


Introduction
Recently, significant research activity has been focused on the so-called core-softened models in which the repulsive part of the intermolecular interaction potential exhibits a softening region at "intermediate" interparticle separations, in addition to usual hard or soft (as in the case of Lennard-Jones (12,6) potential) repulsive branch at short distances. The softening region can be described by a linear or nonlinear ramp, a shoulder, a single or multiple attractive well, or a combination of all these features [1][2][3][4][5][6][7][8][9][10][11]. Model fluids with spherically symmetric core-softened potentials exhibit anomalous thermodynamic and dynamic behaviours [12][13][14] that usually occur in real fluids with directional interparticle interactions, e.g. water, silica, some liquid metals and phosphorus. Examples of the anomalous behavior include the existence of a density maximum as a function of temperature, an increase of the diffusion coefficient upon compression, and, for some systems, the existence of multiple fluid-fluid transitions.
Among different models of the core-softened potentials, the potential used by Barros de Oliveira et al. [22][23][24] that yields both the liquid-gas phase transition and thermodynamic, structural and dynamic anomalies seems to be interesting. The potential is the sum of a Lennard-Jones potential of the well depth ε f and a Gaussian barrier centered on radius r = r 0 with the height of aε f , (1) Depending on the choice of parameters, equation (1) represents the whole family of two length scales intermolecular interactions, from a deep double-well potential [3,25,26] to a repulsive shoulder [4,27,28]. Moreover, due to its continuity, the potential (1) is convenient both for the Monte Carlo and molecular dynamics simulations. For specific choices of the parameters, the potential (1) can exhibit a double well behavior, similar to the potential studied by Cho et al. [3,25]. However, the attractive double well may bring both the liquid-gas phase transition and the anomalies to higher temperatures into an unstable region of the phase diagram in the pressure-temperature plane [26]. In order to circumvent this behavior, several studies [22][23][24][29][30][31][32] were carried out for a specific set of parameters that give rise to a potential with a repulsive ramp followed by a small attractive well. In particular, the following values were often used: a = 5, c = 1 and r 0 = 0.7σ f . It was shown [22][23][24] that for the above set of parameters, the pressure-temperature curves have a minimum if calculated at certain densities (i.e. (∂P/∂T ) ρ = 0). Consequently, the derivative (∂ρ/∂T ) P = 0 for some thermodynamic states. The region of density anomaly corresponds to state points for which (∂ρ/∂T ) P > 0 and is bounded by the locus of points for which the thermal expansion coefficient is equal to zero.
In our recent work [31] we applied grand canonical Monte Carlo simulation and integral equations with hypernetted chain (HNC), as well as Rogers-Young closures to study the above defined model. We found still another anomaly in the system, namely we demonstrated that for some temperatures the derivative of the density with respect to the chemical potential exhibits a minimum at a certain density, followed by a maximum. Also, a peculiarity in the dependence of the specific heat upon density was found. We compared the pair distribution functions resulting from integral equations and from computer simulations and established that common HNC approximation was not successful in capturing the region of anomalies in contrast to Rogers-Young approximation that imposes thermodynamic self-consistency by construction. However, the HNC was very accurate at high fluid densities. It is worth mentioning that a popular mean spherical approximation has not been used to study the model in question because this approach intrinsically requires an adequate splitting of the potential into a short-and long-ranged parts. This approximation does not seem to be beneficial for some soft-core models with two length scales and of this study in particular.
The results of integral equations reported in the literature for the potential (1) were based on the singlet OZ equation [33]. A more sophisticated approach, but one that is also more demanding of computational resources, is based on the inhomogeneous OZ equation [34].
The inhomogeneous OZ equation has usually been used not only to calculate the structure and thermodynamic properties of fluids in contact with surfaces [35][36][37][38][39], but also to determine the structure and the surface tension at the liquid-vapor interface [40]. However, one can also assume that for a single-component fluid, the source of an external potential field is just a single distinguished particle identical to all remaining molecules of the system. This is the so-called Percus' trick [41,42] or the "test particle" method. The method was extended later by Attard [43,44] in the framework of the second-order OZ approach for bulk fluids, see e.g. [45][46][47] for the discussion of related important issues. Moreover, this approach was successfully applied to several model fluids with spherically symmetric intermolecular potentials [48][49][50][51][52][53][54]. Better description of the desired properties or in some cases even qualitatively new findings in comparison to the singlet level theory were obtained. However, to our best knowledge, this approach has never been tested for the core-softened models of the interparticle potentials. Therefore, the motivation of this work is to use the second-order integral equations to the fluid of particles interacting via the core-softened potential (1). The second-order equations are applied to the study of the microscopic structure and thermodynamic properties of the fluid, in particular, in view of its anomalous behavior in a certain region of thermodynamic states. The results are compared with computer simulation data and with the predictions of the singlet integral equations [31]. We explore here the second-order Percus-Yevick and hypernetted closures to the nonuniform OZ equation.

Theory
The OZ equation for an inhomogeneous fluid, wherein the density ρ(r) is not constant reads where h 2 (r 1 , r 2 ) and c 2 (r 1 , r 2 ) are the total and direct correlation functions of an nonuniform fluid. We have used here the subscript "2" in order to distinguish these functions from the common uniform fluid correlation functions. Generally, equation (2) applies when the inhomogeneity is due to an external potential field. However, one distinguished molecule can be also considered as a source of the external potential. In such cases equation (2) can be solved using any of the common approximations, as HNC or Percus-Yevick (PY) approximation In the above u(|r 1 −r 2 |) is the pair potential and β = 1/kT is the inverse temperature. However, we stress that the above closures do not relate uniform, but rather nonuniform correlation functions and therefore we refer to the results of equation (2) with the closure given by equation (3), or by equation (4), as to the HNC2 and PY2 results, respectively, in contrast to the usual HNC and PY theories [33].
To have a set of equations completed we also need a relation between ρ(r) and the pair correlation functions. The exact equation developed by Lovett, Mou, Buff and Wertheim (LMBW) [55,56] reads where y(r 1 ) is the one-particle background (cavity) function, ρ(r 1 ) = exp[−βv(r 1 )]y(r 1 ) and v(r 1 ) is the potential due to the particle which is regarded as the source of the inhomogeneity. In our case, v(r 1 ) is just the pair potential, v(r 1 ) ≡ u(|r 1 − 0|), where the distinguished particle is set at the origin, r 2 = 0. One should mention here an important difference between the theories used to study systems in contact with a surface [35][36][37][38][39] (or the systems involving a gas-liquid interface [40]) and theories based on the method of Attard. In the former case the two-particle correlation functions reduce to the functions for bulk fluids, providing the both particles to be located far away from an inhomogeneous region. In the Attard's approach, the function h 2 corresponds rather to the thirdorder (conditional) correlation function of a bulk fluid, while the "usual" pair correlation function is related to the local density. Indeed, because the bulk pair distribution function, g(r), gives the probability density of finding a pair of fluid molecules at a separation, r, the g(r) is related to the local density, ρ(r), via where |r| is the distance from the distinguished particles that is the source of the system inhomogeneity and ρ = lim |r|→∞ ρ(|r|) is the bulk density. Of course, the cavity function y(r) that enters equation (5) satisfies the relation g(r) = exp[−βu(r)]y(r).
To solve equations (2) and (5) with the closure given by equation (3) or equation (4) we use the numerical algorithm Attard's [43] that relies on the expansion of the two-particle functions (h 2 and c 2 ) in series of Legendre polynomials, for details see [43].
Note that the closure equations (3) and (4) constitute an approximation, similar to the common bulk theory, whereas the OZ equation (2) and the equation for the profile, (5) are exact.
The knowledge of the total, h(r), and the cavity, y(r), functions enable us to determines the bridge function, B(r). The bridge function is defined as, see e.g. [33,57] where γ(r) = h(r) − c(r). The bridge function plays a key role in theories based on the singlet OZ equation (7). It is given by the smallest set of diagrams in the graphical expansion, and all other functions can be calculated from it via the OZ equation [33].
The main difference between the first-and the second-order integral equation theories for bulk fluids is that in the latter approach, the bridge function results from the calculations performed, whereas in the case of first-order theories the bridge function is imposed as a closure. Previous calculations carried out for hard spheres, as well as for Lennard-Jones fluids, gave insight into the course of the bridge functions. However, to our best knowledge no results of calculations of the bridge function for core-softened potential models have been presented so far.

Results and discussion
The second-order integral equations were solved using the Attard's algorithm [43]. Usually, 80 Legendre polynomials were used and the grid in r was 0.04σ. To test the accuracy we increased the number of Legendre polynomials to 120 and decreased the grid size to 0.025σ. Moreover, we also solved the singlet integral equations and carried out grand canonical ensemble Monte Carlo simulations. The details on the two latter methods are presented in our previous work [31]. We also note that the pressures given below were calculated from the virial equation of state (see equation (6) of [23]).
Similarly to our previous work [31], all the calculations have been carried out assuming that the parameters of the potential (1) are a = 5, c = 1 and r 0 = 0.7σ f . First, we concentrate on a comparison of the pair distribution functions obtained from the second-order OZ equation, the singlet OZ equation and the Grand Canonical ensemble (GCMC) computer simulations [31]. Comparisons are for thermodynamic conditions that cover the density anomaly region found in molecular dynamics studies [23]. In figures 1 (a)-(c) we display the functions g(r) evaluated at the density ρ * = ρσ 3 f = 0.1 and at three temperatures, T * = kT /ε f = 0.5, 0.3 and 0.2. At the highest considered temperature the results of the singlet theories (PY, HNC and Rogers-Young, (RY)), the second-order theories, PY2 and HNC2 are compared with GCMC data. Before discussing the results shown in figure 1, we recall that the singlet RY closure contains one adjustable parameter and this parameter is evaluated requiring the equality of pressures from the virial and compressibility routes, for details see [31].
In figure 1 (a) we see that the singlet PY approximation yields the radial distribution function that significantly deviates from all remaining results. The second-order PY2 approximation is definitely more accurate than the singlet PY theory, but its results still differ much from computer simulations. In contrast to the singlet PY theory, the singlet HNC closure yields the results that qualitatively agree with computer simulations, but among all the singlet theories the thermodynamically self-consistent RY closure gives the best predictions. For the thermodynamic state from figure 1 (a) the second-order hypernetted chain approximation reproduces the GCMC data with the highest precision. We stress that in contrast to the singlet RY theory the HNC2 approach contains no adjustable parameter.
The accuracy of the second-order HNC2 theory is also confirmed by the results presented in figure 1 (b). HNC2 theory well reproduces the location and the height of the first maximum of g(r) and the shoulder at distances r/σ f < 1.4. However, the first minimum of g(r) at r/σ f ≈ 3.32 is a bit better predicted by the RY approximation. Similar conclusions can be drawn from the inspection of the curves presented in figure 1 (c). Since the PY2 approximation actually fails to reproduce the simulation data with a reasonable accuracy, the PY and PY2 results are omitted in figure 1 (c) and in the following figures 2 and 3. The singlet PY approximation is quite accurate when applied to a hard-sphere system [33]. On the contrary, in the case of soft potentials, e.g.,

13601-4
Core-softened model fluid Gaussian-like potential, the singlet HNC approximation works quite satisfactory [58]. It is thus not surprising that for core-softened potential, the hypernetted chain approximation (at singlet and second-order levels) performs better than the the Percus-Yevick closures. However, the observed big differences between the predictions of these approximations are surprising, because in the case of the systems studied so far [48][49][50][51][52][53] the differences between the PY2 and HNC2 approximations were rather small [59]. Further tests of the accuracy of the HNC2 approximation are shown in figures 2 and 3. Figure 2 shows the results at the density ρ * = 0.14, consecutive parts ad are for the temperatures T * = 0.5, 0.3, 0.2 and 0.15, respectively. Again, the HNC2 approximation is superior over the best singlet RY theory. In particular, it better describes the formation of the first peak of g(r) within the repulsive ramp of the potential (1) (the plot of the potential (1) gives figure 1 of [23]), though both RY and HNC2 theories underestimate the height of this maximum. By contrast, the singlet HNC approach describes much higher first peak of g(r) compared to GCMC data.
Finally, we calculated the radial distribution functions at higher densities, ρ * = 0.16 and 0.20 ( figure 3). For these densities, the second-order HNC2 equation performs well. In particular, its predictions of the height of consecutive maxima and minima of g(r) are more accurate than the singlet RY approximation. However, for some state points, the accuracy of the singlet HNC equation  is quite good and even better than the HNC2 theory, cf. figures 3 (b) and 3 (c). Before discussing an anomalous behavior of the system in the pressure-temperature plane, that results from the second-order integral equations, we briefly recall previous findings. The line of temperatures of maximum density (TMD) was determined from NPT molecular dynamics [23] and singlet PY [23] and RY [23,31] integral equations. It was also demonstrated [23,31] that the singlet HNC approach fails to predict this anomaly.
Molecular dynamic simulations revealed that for densities 0.12 < ρ * < 0.14, the pressuretemperature curves at constant density have a minimum (see also our comment in the introductory section) which implies density anomaly. The maximum temperature at which the density anomaly still exists is T * max ≈ 0.25 [23]. The predictions of the singlet RY approximation depend on the method according to which its adjustable parameter was determined [31]. Barros de Oliveira et al. [23] obtained the adjustable parameter of the RY equation by checking the consistency between the compressibilities, calculated from the virial and compressibility equations of state. However, in our work [31] we imposed the so-called global coexistence criterion, according to which the pressures from two above mentioned thermodynamic routes were compared. Barros de Oliveira et al. [23] found that the density anomaly exists for the densities from the range of [0.12, 0.14] and that T * max ≈ 0.23. Our calculations, however, led to a wider interval of the densities, [0.12, 0.19]. We should also note that the singlet PY approximation predicts [23] much wider region of densities, at 13601-6 Core-softened model fluid 2 4 r/σ f which the density anomaly occurs, [0.13, 0.3]. Instantaneously, this region is shifted towards much higher temperatures, with T * max = 0.86. Figure 4 shows the isochores evaluated from different theories (the reduced pressure is defined as P * = P σ 3 f ε f ). At low densities (ρ * < 0.07; the relevant curves are not shown for the sake of brevity) none of the approximations predicts the existence of the density anomaly. For ρ * = 0.08 (part a) only the PY2 theory yields a weak minimum of the pressure, all the remaining approximations do not show any anomaly. However, it is interesting to note that the results of the singlet PY and HNC2 approximations almost coincide, except for low temperatures. For ρ * = 0.1 (part b) only the PY2 isochore exhibits anomaly, but now the minimum of the pressure is preceded by a maximum. This behavior is completely unexpected, but, it is rather an artifact of the theory, because no computer simulations [23] provide its confirmation. For ρ * = 0.14 (part c) all the approximations, but the HNC1 predict the density anomaly, whereas at higher bulk density, ρ * = 0.2, (part d) only PY1 and PY2 predict anomalies.
The range of the densities for which the HNC2 leads to anomalous behavior is rather narrow, 0.13 ρ * 0.155. As we have shown in our previous work [31], the RY approximation with the global coexistence criterion predicts anomalous behavior up to ρ * ≈ 0. 19. Surprisingly, the PY2 approximation yields two density anomaly regions: the first one is for the densities [0.07, 0.11] and the second one is for ρ * 0.13. We must recall, however, that the structure predicted from PY2 13601-7 approximation differs very much from the computer simulation results. Summary of our calculations is presented in figure 5. We show here TMD lines resulting from different approximate theories.
Since the PY2 predictions quantitatively differ from the computer simulation results, the relevant curve has been omitted. The TMD line from the RY approach evaluated using global coexistence criterion differs from that obtained by Barros de Oliveira et al. [23] None of the approximate theories is capable of reproducing the simulation results at quantitative level. Finally, we show the results of the bridge function calculations. The bridge function for the model in question has never been investigated so far. The function γ(r), entering equation (8) was calculated from whereh(k) is the Fourier transform of the function h(r). We should stress that the calculations of the Fourier transformh(k) and the inverse transform in equation (9) should be carried out with a special care, as described by Kolafa et al. [60]. Our calculations were performed at a constant density equal to ρ * = 0.152 for a set of temperatures across the TMD line, as well as at a constant temperature, T * = 0.175 and for several densities, again across the TMD line (cf. figure 5). The functions B(r) significantly differ from the functions for hard-spheres [48,60]. In particular, the decay of the bridge functions for the core-softened potential is much slower than for hard-spheres. Even for r/σ f ≈ 10, the bridge functions exhibit well pronounced oscillations. These oscillations almost vanish at a distance as large as r/σ f ≈ 16, whereas for hard-spheres at very high densities, the oscillations already vanish at r/σ f ≈ 5, cf. figure 3 of [60]. Crossing the TMD line seems to have no effect on the shape of the bridge functions. The evolution of the functions B(r) along the temperature and the density branches ( figure 6) is smooth and we do not observe any peculiarities connected with specific thermodynamic behavior of the system. However, the range of oscillations of a bridge function is worth stressing once again. Let us briefly summarize our findings. We have demonstrated that the radial distribution functions predicted by the second-order HNC2 approach are more accurate than the predictions of any singlet theory. In contrast to the singlet RY approach, the HNC2 approximation does not involve any adjustable parameter. However, the PY2 approximation is inaccurate and the distribution functions resulting from it significantly differ from those of all remaining theories. The discrepancies between PY2 and HNC2 approximations are puzzling and we cannot offer any physical interpretation of these trends. The second-order HNC2 approach yields the density anomaly in the system  for densities from the interval [0.13, 0.15]. The maximum temperature at which this phenomenon is observed is T * max ≈ 0.2. Moreover, for the first time we have evaluated the bridge functions for the core softened fluid. Unlike the hard-sphere or Lennard-Jones fluids, the oscillations of the bridge functions extend over much larger interparticle separations. A successful parametrization of the bridge function similar to hard-sphere systems, see e.g. [60][61][62], would be helpful in developing new closures for singlet level theory for a class of core-softened potentials.