Density-dependent interactions and thermodynamic consistency in integral equation theories

In this paper we present an alternative formulation of the well-known integral equation approximations designed to keep a consistent approach to the determination of thermodynamic properties in the case of density-dependent interactions. Obviously, residual inconsistencies inherent to the approximate character of the closure relations of the Ornstein–Zernike equation will not be corrected. In this connection, we will show how this approach is particularly successful when applied in conjunction with approximations in which the aforementioned inconsistencies are minimal, as is the case of the optimised Reference Hypernetted Chain equation. As a case study we will consider the Derjaguin–Landau–Verwey–Overbeek model of charged colloids which is one of the simplest realisations of density-dependent interactions.


Introduction
Density-dependent potentials are ubiquitous in the field of liquid state and soft matter physics. Such is the case of colloidal dispersions [1], in which the most widespread interaction model is described by the theory of Derjaguin-Landau-Verwey-Overbeek (DLVO) [2,3]. In this particular instance large charged colloidal particles are immersed in an electrolyte solution, experiencing a screened Coulomb interaction that depends on the density of the counterions (and consequently on the density of the colloidal particles). Also, one of the most representative cases of densitydependent interactions is that of the effective ion-ion potential in liquid metals when treated in the nearly free electron model [4]. Here, the density-dependence appears in the ion-ion interaction through Fermi's wave number, k F . This quantity is a simple function of the electron density, which once again is directly related to the ionic density through the electroneutrality condition. Other effective interactions in metals can be found in the embedded atom method [5,6], the Finnis-Sinclair potentials [7], or the effective medium theory [8], among many others. In all cases, effective potentials can be thought of as the result of a coarsegraining procedure, in which a certain number of degrees of freedom are integrated out to yield a much simpler pairwise additive potential. This coarsegraining can be done explicitly when deriving an interaction potential in a theoretical treatment (as is the case of the DLVO or the effective ion-ion interaction in metals). Additionally, effective densitydependent potentials can be generated from inversion procedures applied to scattering or ab initio simulation data [9][10][11]. In this case, for each thermodynamic state a different potential will be generated, and the coarse-graining of various degrees of freedom is implicitly performed.
Useful as they are, these types of effective interactions must be dealt with some care. To the best of our knowledge, Ascarelli and Harrison [12] were the first to point out that the virial expression must be modified to incorporate explicitly the effect of the density dependence of the potential on the evaluation of the pressure. More recently [13], it was stressed that the usual form of the fluctuation theorem is no longer valid when dealing with effective density-dependent potentials, being inconsistent with the virial theorem result. Besides, simulation procedures for this type of interaction must also be adapted when sampling inhomogeneous systems and/or systems in the vicinity of first order phase transitions [14]. In a recent work by Louis [15], it was argued that the inconsistency between virial and fluctuation theorem routes stems from the fact that the virial route, defined in the canonical ensemble, samples a single density, whereas the usual derivation of the fluctuation theorem isothermal compressibility is made in the grand canonical ensemble. Therefore, in the latter instance all densities are sampled, and this implies that different values of the interaction are also sampled. This precise feature is disregarded when the standard expression for the isothermal compressibility is used. Additionally, in [15], it was also stressed the approximate character of the coarse-graining procedures. This implies that depending on the nature of the original problem that is being reduced to a system of 'effective particles' interacting via pairwise additive potentials and the route followed to perform the reduction, one will get different expressions for the same quantity. These different paths to calculate thermodynamic quantities will agree more or less depending on how severe the coarse-graining approximations have been for a given system. Now, the behaviour of the effective systems mentioned above, can be analysed in terms of computer simulation -if needed, using local density approaches [14] -perturbation theory or integral equation methods. In the latter instance, one of the most powerful approaches is the implementation of thermodynamically self-consistent closures [16,17]. A well known drawback of the integral equation approximations is the lack of consistency in the thermodynamic properties derived from different routes (energy, virial, free energy or fluctuation theorem compressibility routes). This lack of consistency stems from the approximate character of the closure relation imposed on the Ornstein-Zernike (OZ) equation. In this context, the Hypernetted Chain (HNC) equation is known to be fully consistent except for the virial-fluctuation theorem compressibilities [18]. This residual inconsistency can be further minimised incorporating a reference system approximation to the bridge function, which is neglected in the HNC approach. This constitutes the Reference Hypernetted Chain (RHNC) approximation [19]. In the case of density-dependent interactions, the consistency between the energy, virial and free energy routes is well preserved if the correct virial expression [12] is used. Now, the self-consistent approaches exploited in [16,17] are based on the use of parameterised hybrid closures, in which the closure parameters are tuned in order to enforce consistency between virial and fluctuation theorem compressibilities. Taking into account the discussion of the previous paragraph it is obvious that in principle this would not be a feasible task in the case of density-dependent pair potentials, since the expressions for the two quantities we intend to make consistent are derived on a different footing [15]. The aim of this paper is to present a reformulation of the integral equation approaches that, from an operational standpoint, make possible the definition of consistent virial pressures and fluctuation theorem compressibilities for density-dependent interactions. These quantities should be fully consistent whenever an exact closure relation to the Ornstein-Zernike equation were available.
The rest of the paper can be sketched as follows. In the next section we review the essentials of the HNC and RHNC integral equations and reformulate a corresponding set of integral equations that should lead to consistent thermodynamics. Note that here, when referring to consistent thermodynamics we do not take into account the residual inconsistency introduced by the approximate closure relations in use, i.e. the inconsistency we intend to cure is the one derived from the use of density-dependent potentials in the definition of quantities such as virial pressure and fluctuation theorem isothermal compressibilities. A new integral equation for density-dependent interactions is thus formulated, and general expressions to determine pressures and isothermal compressibilities are derived. Then in Section 3, we apply the formalism to the DLVO potential and present our most significant results and conclusions.

Effective interactions and consistent integral equations
A fundamental approach of the theory of simple fluids relies on the use of integral equations. In this type of approach, the pair distribution function g 0 (r; , T ), with r the interparticle distance, the average number density and T the temperature, is calculated from a given intermolecular pair potential V 0 (r). In the specific case of OZ-type methods, one starts with the exact OZ relation relating the total correlation function h 0 (r; , T ) ¼ g 0 (r; , T ) -1 and the direct correlation function c 0 (r; , T ). An approximate equation (a closure relation) which links both functions through the intermolecular pair potential is needed. In general, this closure can be written as where ¼ 1/k B T with k B the Boltzmann constant, and b 0 (r; , T ) is the so-called bridge function. One of the simplest approaches implies the neglect of this term (b 0 ¼ 0), and this renders the HNC equation.
The RHNC goes a step beyond, by approximating the unknown bridge function by one of the well-defined reference systems . The most widely used reference system in the case of simple fluids is obviously the hard-sphere fluid, for which the bridge function can easily be calculated [20,21]. Once the non-linear integral equation resulting from Equations (1) and (2) is solved, the thermodynamics can be derived from the Helmholtz freeenergy per particle f 0 (, T ) which can be expressed in two exact and equivalent ways. In the pair correlation (or g-) route the thermodynamic potential is given by where Ã is the thermal de Broglie wavelength and g 0 (r; , Tj) is the pair correlation function for the fluid with pair potential V 0 (r), whereas in the direct correlation (or c-) route reads The main difference between these two routes is that in Equation (3) the pair correlation function is charged with interactions (at constant density) and V 0 (r) appears explicitly, whereas in Equation (4) the direct correlation function is charged with density (at constant interactions) and there is no explicit dependence on the pair potential at all. Although the direct correlation route is independent of the particular form of the interactions, which may not be pairwise additive, in order to compare both routes, the pair potential V 0 (r) has to be introduced in the closure relation as in Equation (2). Let us now consider a fluid interacting with a density-dependent pair potential V(r; ). At first sight it should be concluded that all we have to do in this case is to replace in Equation (3) V 0 (r) by V(r; ), while keeping Equation (4) unchanged. We have nevertheless to be cautious with this naive generalisation. Note that when we replace V 0 (r) by V(r; ) all the structural functions, and consequently, all the thermodynamic quantities, acquire an additional dependence on the density. When this additional density dependence (indicated in what follows by a subscript ) is taken into account in the pair correlation function, say g (r; , T ), Equation (3) can be easily generalised yielding: where g (r; , Tj) is the pair correlation function for the fluid with pair potential V(r; ). Within this route, the standard integral equation methods given by Equations (1) and (2) can be extended with the substitutions V 0 (r) ! V(r; ), g 0 (r; , T ) ! g (r; , T ), and c 0 (r; , T ) ! c (r; , T ), where c (r; , T ) is the direct correlation function for the fluid with potential V(r; ).
A straightforward transformation of Equation (5) within the HNC leads to: with and wherec ðk; , T Þ andh ðk; , T Þ are the Fourier transforms of c (r; , T ) and h (r; , T ), respectively. Equations (6)(7)(8) involve only the average number density of the fluid. A similar expression can be derived in the RHNC approximation. When the hardsphere fluid is taken as the reference system [19] the derivation of Equations (6-8) for the free energy in the RHNC is straightforward.

A consistent approach
When the integral equation methods are applied using the direct correlation route we are faced with an unusual feature. At first sight, one could be tempted to generalise Equation (4) as: since the c-route seems not to depend on the pair potential. A cumbersome but straightforward calculation shows however that for an analytic densitydependent pair potential, the density expansions of Equations (5) and (9) differ even at the level of the third virial coefficient! This can be shown by first considering the (formal) virial expansions of g (r; , Tj) and c (r; , T ), whose virial coefficients can be further expressed in terms of convolutions of Mayer functions evaluated with V(r, ) and V(r, ), respectively. To obtain the full density expansions of the excess free-energy, one still has to introduce the virial expansions of these density-dependent virial coefficients. This thermodynamic inconsistency can be explained as follows.
Note that the cornerstone of the theory of simple fluids is the equivalence of Equations (3) and (4). In the latter the parameter appearing in c 0 (r; , T ) operates on the density at constant interactions. So, in order to maintain the equivalence of the two routes for density-dependent pair potentials, the density appearing in V(r; ) has to be considered as a constant (i.e. constant interactions) when the direct correlation function is charged with density. This would remove the inconsistency that stems from sampling a single global density in the case of Equation (3) and sampling all densities in Equation (4). This lack of consistency is identical to the one highlighted in [15] between virial pressure and fluctuation theorem compressibilities. Now, in other words, we can generalise Equation (4) as: which with the change of variable 0 ¼ , reads In this case, we have to solve the Ornstein-Zernike equation together with the HNC/RHNC closure Notice the -dependence (and not 0 ) of V(r; ).
In summary, for a prescribed average number density , we thus obtain a non-linear integral equation involving two different densities and 0 with 0 0 . The variable 0 has to be discretised and the set of integral equations has to be solved for a series of 0 i values, yielding a set of c (r; 0 i , T ) functions that can be numerically integrated to yield the Helmholtz free-energy through the direct correlation route.
Observe that in the original (inconsistent) formulation, the c-route to the Helmholtz free energy reads where now c 0 (r; 0 , T ) will be evaluated for V(r; 0 ) as well.

Effects on the virial pressure and the isothermal compressibility
Now, if we take into account that where P is the pressure, taking the density derivative of Equation (11), one gets The first term of the right-hand side of Equation (16) is precisely the integral of the inverse isothermal compressibility as calculated from the fluctuation theorem in the case of density independent potentials. Now by comparing Equation (16) with the virial equation as derived by Ascarelli and Harrison [12] P ¼ À 1 2 2 Z dr g ðr; , T Þ r 3 @Vðr; Þ @r where the first two terms of the right-hand side of Equation (17) are retained in the case of density independent potentials, the consistency of Equation (16) and Equation (17) -subject to the unavoidable inaccuracies of the closure relation -implies that i.e. the terms P that explicitly incorporate the effect of the density dependence of the interaction on the pressure should be identical when calculated from the c-and g-routes. Now the question would be how to determine the density derivative of the direct correlation function. The simplest way to do this is to follow Belloni's method [22] and differentiate both the OZ and the corresponding closure relation with respect to the density (not 0 ). This leads to a rapidly convergent linear integral equation that can be solved immediately after the HNC solution has been numerically determined. If one proceeds further to differentiate Equation (16) with respect to the density, one gets [23,24] in which, again, the first term of the right-hand side is the isothermal compressibility of a system interacting via density independent potentials. The remaining terms should make this expression consistent with the isothermal compressibility obtained by differentiation of the virial equation, Equation (17).

A case study: the DLVO potential
As a test case, we will apply the above integral equations to the DLVO potential. This reads where V HS (r/) is a hard-sphere potential of diameter , -Ze the charge, the dielectric constant and 2 ¼ 4 B ( B ¼ e 2 / B T ). With the definition of x ¼ r/ and R being the radius of the colloidal particles, the interaction can be expressed as Here we have defined The calculations presented below have been carried out for b ¼ 1, R* ¼ 1, and T* ¼ B T/" ¼ 4. The integral equations have been solved on a discretised grid with Ár ¼ 0.01 and 4096 points. In the RHNC calculations the hard-sphere reference diameter has been optimised following the prescription of Lado et al. [19], requiring the minimisation of the free energy.
In Figure 1 we depict the density dependence of the free energy, calculated in the HNC using the g-route -Equation (6) -and the c-route -Equation (14). The lack of consistency of both results is evident, even for densities as low as 3 ¼ 0.1. Now, when the consistent c-route -Equations (11-13) -is used, the situation is substantially different. The agreement when comparing the results using both routes is similar to the one found when dealing with density independent potentials. The residual lack of consistency can be attributed to the approximate character of the closure. This fact can be further confirmed by a careful analysis of Figure 2, in which the same quantities, calculated in the RHNC approximation, are plotted. It is known that when considering the lack of consistency between the virial and fluctuation theorem compressibilities, the RHNC somewhat improves upon the HNC [19]. This improvement translates in our case to a better agreement between the g-and the consistent c-route -Equation (11) -to thermodynamics, and this is reflected in Figure 2. One readily appreciates that the consistent c-route results (filled circles) fall on top of the g-route curve (solid line).
We can now analyse to what extent the previous results affect the calculation of the pressure. This is best illustrated comparing the quantities which are plotted in Figure 3 for the HNC approximation. Calculations in the RHNC for Equation (22) are somewhat more involved since they require the density derivative of the reference system bridge function in order to calculate the derivative of the pair  Figure 1. Free energy calculated in the HNC approximation via the g-route (solid line), the c-route (dotted line) and the consistent Equations (11-13) (filled circles) for the DLVO potential with the parameters indicated below Equation (21) in the text. The agreement between the g-route and the consistent equation is notorious. correlation function. This is numerically cumbersome and therefore for our purposes we have limited ourselves to the HNC level of approximation. We see in Figure 3 that the agreement between the new c-route density-dependent term with the term derived by Ascarelli and Harrison [12] is fairly good up to moderate densities. Deviations at higher densities, once more, are due to the approximations implicit in the HNC closure. As was seen in the case of the free energy calculation, the use of the RHNC in this instance is expected to reduce the residual inconsistency.
In summary, we have shown how two well-known integral equations can be reformulated so as to render the c-and g-routes to thermodynamics consistent for density-dependent potentials. This consistency, which is limited by the approximations implicit in the closure relation, can be improved thanks to the implementation of more elaborate closures. Thus, in principle, it should be possible to build hybrid closures along the lines of [16,17] and enforce the virial-fluctuation theorem compressibility consistency. Given the form of Equations (11)(12)(13), a global consistency approach [25] -i.e. using the pressures and not the compressibilities -should be easier to implement.  Figure 2. Free energy calculated in the RHNC approximation via the g-route (solid line), the c-route (dotted line) and the consistent Equations (11-13) (filled circles) for the DLVO potential with the parameters indicated below Equation (21) in the text. The agreement between the g-route and the consistent equation has improved with respect to Figure 1.  Figure 3. Evolution of the terms that account for the density dependence of interaction in the evaluation of the pressure (see Equations (16)(17)), as calculated from the g-route (Equation 23), and c-route (Equation (22)). These calculations have been carried out in the HNC approximation for the DLVO potential with the parameters indicated below Equation (21) in the text.