text:
# Different thermodynamic pathways to the solvation free energy of a spherical cavity in a hard sphere fluid

## I Introduction

The solvation free energy determines how readily a solute can be dissolved in a given solvent fluid. This plays an important role in many chemically and biologically important processes, perhaps most notably in hydrophobic interactions in water. A significant part of the solvation free energy arises from the required expulsion of solvent molecules from the region occupied by the harshly repulsive molecular core of the solute. These very strong “excluded volume” interactions can significantly perturb the local density around the solute and cause simple approaches based on gradient expansions to fail.

These effects can be seen most clearly in the simple model system treated in this paper. We will calculate the excess or solvation free energy associated with the insertion of a spherical cavity with radius $`R_v`$ into a hard sphere fluid, whose molecules have diameter $`\sigma `$. By definition, the centers of the solvent molecules are completely excluded from the region of the cavity, which thus acts like a hard core external field.

This system has many interesting limits. When the exclusion field or cavity radius $`R_v`$ equals $`\sigma `$, then the cavity acts like another solvent particle and the solvation free energy is directly related to the chemical potential of the solvent. As the cavity radius tends to infinity it effectively turns into a hard wall and the relevant thermodynamic quantity is the surface free energy or surface tension associated with a hard wall in a hard sphere fluid. A cavity with radius $`R_v=\sigma /2`$ acts like a hard core “point solute” of zero diameter. Even shorter-ranged hard core fields or “tiny cavities” with $`R_v\sigma /2`$ are also of interest, since the induced structure and solvation free energy of a tiny cavity can be calculated exactly. This limit can thus serve as a nontrivial check on approximate methods.

The most commonly used method today for such problems is weighted density functional theory (DFT) evans:density . Here one attempts to describe the free energy directly as a functional of some kind of smoothed or weighted average of the nonuniform and often rapidly varying singlet density. This has the advantage that the free energy is obtained directly and by construction the associated fluid structure (obtained by functionally differentiating the free energy) is consistent with the approximate free energy. However the choice of appropriate weighting functions is by no means obvious and a number of different and often highly formal schemes have been proposed.

We focus instead in this paper on making direct use of structural information about the nonuniform solvent fluid to obtain the solvation free energy. We believe this allows physical intuition to play a more central role and we can take advantage of the recent development of a generally very accurate theory relating the structure of a nonuniform hard sphere fluid to the associated external field hlr .

As we will see below, the free energy can then be calculated by integration, starting from an initially known state (e.g., the uniform fluid) and determining the free energy changes as the solute-solvent interaction (the hard core external field) is “turned on”, or alternatively, as the density is changed from the initial to the final state. There exist many possible routes from the initial to the final state, and we will generally refer to them as thermodynamic pathways. If exact results are used for the intermediate values of the structure and associated fields, then all these different pathways will give the same (exact) result for the free energy.

In practice, of course, approximations will have to be made and the different pathways will generally yield different results. This is sometimes referred to as the “thermodynamic inconsistency” of structurally based methods evans:density . But this can be viewed more positively as giving one the freedom to choose particular pathways that could be relatively insensitive to the errors that exist in the structural theory, and we will try to use this flexibility to obtain the most accurate results. Moreover, there is an inherent smoothing of the structural information in the integration used to obtain the free energy. The differences in free energy predicted by different pathways will also give us some indication about the overall quality of the theory.

This approach generally requires the density profiles and associated fields of all the intermediate states along the various pathways, and thus a fast and accurate method for determining these quantities is crucial for computation efficiency. We will use here the generally accurate hydrostatic linear response (HLR) equation hlr  proposed by Katsov and Weeks. A different physically motivated derivation of the HLR equation is given below.

We will also examine the alternative free energy predictions that arise from a theory closely related to the HLR equation, the Gaussian field model (GFM) developed by Chandler gft . For a solute with a hard core the GFM proposes an approximate partition function from which the associated density response can be derived. In the particular case where a rigid cavity is inserted into a hard sphere fluid, the HLR and the GFM approaches turn out to make identical predictions for the induced structure. Thus structurally based routes to the free energy involving only hard core fields will give the same results. In addition, one can use the approximate GFM partition function to evaluate the solvation free energy directly. However, as we will show later, the latter approach tends to produce less accurate results. This deficiency shows up even more strongly in the tiny cavity limit where the structural predictions of the HLR and the GFM are exact, and several pathways giving the exact free energy can be found. This illustrates the advantage of considering a variety of thermodynamic pathways that can make best use of the available structural information.

## II Density response to an external field

### II.1  The HLR equation

We describe the system using a grand canonical ensemble, and thus want to determine the excess grand free energy arising from insertion of a spherical cavity or hard core external field of radius $`R_v`$. To derive the HLR equation hlr  we start with the basic linear response equation hansenmac  for a nonuniform hard sphere system in a general external field $`\varphi (𝐫)`$, with chemical potential $`\mu ^B`$, inverse temperature $`\beta =(k_BT)^1`$ and associated density $`\rho (𝐫;\mu ^B,[\varphi ])`$ $`\rho (𝐫)`$:

$$\beta \delta \varphi (𝐫_1)=𝑑𝐫_2\chi ^1(𝐫_1,𝐫_2;[\rho ])\delta \rho (𝐫_2).$$
(1)
This relates small perturbations in the density and field through the (inverse) linear response function

$$\chi ^1(𝐫_1,𝐫_2;[\rho ])\delta (𝐫_1𝐫_2)/\rho (𝐫_1)c(𝐫_1,𝐫_2;[\rho ]).$$
(2)
Here $`c(𝐫_1,𝐫_2;[\rho ])`$ is the direct correlation function of the nonuniform hard sphere system. The notation $`[\rho ]`$ indicates that these correlation functions are nonlocal functionals of the density $`\rho (𝐫).`$

Since we want to focus on the effects of the perturbing field, we have used the inverse form of linear response theory prattotf  in Eq. (1), where the field appears explicitly only on the left hand side, evaluated at $`𝐫_1.`$ This provides many advantages in dealing with large field perturbations, as will soon become apparent. In most cases we will consider perturbations about a uniform system with chemical potential $`\mu `$ and density $`\rho (\mu )\rho (𝐫;\mu ,[0])`$. When using this simplified notation $`\rho (\mu )`$ should not be confused with $`\rho (𝐫)\rho (𝐫;\mu ^B,[\varphi ])`$. Similarly, we will let $`\mu (\rho )`$ denote the chemical potential of the uniform fluid as a function of density $`\rho `$. In a uniform system the direct correlation function $`c`$ will take the simple form $`c(r_{12};\rho ),`$ where $`r_{12}|𝐫_1𝐫_2|`$.

But how can we use Eq. (1) to describe the density response to a large field perturbation such as the hard core field of interest here? This linear relation between a (possibly infinite) external field perturbation on the left hand side and the finite induced density change on the right must certainly fail for values of $`𝐫_1`$ where the field is very large. Conversely, Eq. (1) should be most accurate for those values of $`𝐫_1`$ where the field is small — in particular where the field vanishes — and then through the integration over all $`𝐫_2`$ it relates density changes in regions where the field vanishes to density changes in the regions where the field is nonzero.

To treat large fields, we note that for any given $`𝐫_1`$ we can locally impose the optimal condition that the field perturbation vanishes by introducing a shifted chemical potential

$$\mu ^{𝐫_1}\mu ^B\varphi (𝐫_1),$$
(3)
and a shifted external field

$$\varphi ^{𝐫_1}(𝐫)\varphi (𝐫)\varphi (𝐫_1).$$
(4)
Since there is an arbitrary zero of energy and a constant external field acts like a shift of the chemical potential in the grand ensemble, we make no physical changes if we shift both functions by the same amount. In particular $`\rho (𝐫;\mu ^B,[\varphi ])=\rho (𝐫;\mu ^{𝐫_1},[\varphi ^{𝐫_1}]).`$

The superscript $`𝐫_1`$ in $`\mu ^{𝐫_1}`$ indicates a particular value of the chemical potential, which from Eq. (3) depends parametrically on $`𝐫_1`$ through the local value of the field. When $`\varphi (𝐫_1)`$ vanishes, then $`\mu ^{𝐫_1}`$ reduces to $`\mu ^B.`$ We define $`\rho ^{𝐫_1},`$ the hydrostatic density, by

$$\rho ^{𝐫_1}\rho (𝐫;\mu ^{𝐫_1},[0])=\rho (\mu ^{𝐫_1}).$$
(5)
Thus $`\rho ^{𝐫_1}`$ is the density of the uniform fluid in zero field at the shifted chemical potential $`\mu ^{𝐫_1}`$; equivalently $`\rho ^{𝐫_1}`$ satisfies

$$\mu (\rho ^{𝐫_1})=\mu ^{𝐫_1}=\mu ^B\varphi (𝐫_1).$$
(6)
Thus far, we have have merely introduced an equivalent (and apparently more complicated!) way of describing the system in terms of a shifted field and a shifted chemical potential. However this perspective immediately suggests a very simple first approximation to the density response to a slowly varying external field. Since $`\varphi ^{𝐫_1}(𝐫)`$ by construction vanishes for $`𝐫=𝐫_1`$, if $`\varphi ^{𝐫_1}(𝐫)`$ is sufficiently slowly varying, then the region around $`𝐫_1`$ within a correlation length is essentially in zero field. In that case the uniform hydrostatic density $`\rho ^{𝐫_1}`$ is clearly a good approximation to $`\rho (𝐫_1)`$, the exact induced density at $`𝐫_1`$. Moreover, when the field is more rapidly varying, it is natural to introduce a second and even more accurate approximation to the density response.

The hydrostatic density $`\rho ^{𝐫_1}`$ takes account only of the local value of the field at $`𝐫_1`$ by a shift of the chemical potential. The HLR equation improves on this “local field” approximation by using linear response theory to determine the density change from the hydrostatic density induced by nonlocal values of the shifted field $`\varphi ^{𝐫_1}(𝐫)`$. Thus starting from the uniform density $`\rho ^{𝐫_1}`$, we assume a linear response to the shifted field, replacing $`\chi ^1(𝐫_1,𝐫_2;[\rho ])`$ by $`\chi ^1(r_{12};\rho ^{𝐫_1})`$ in Eq. (1) and setting $`\delta \varphi (𝐫)=\varphi ^{𝐫_1}(𝐫)`$ and $`\delta \rho (𝐫_2)=\rho (𝐫_2)\rho ^{𝐫_1}`$. Then the left side of Eq. (1) vanishes (giving the optimal linear response condition), and we have

$$0=𝑑𝐫_2\chi ^1(r_{12};\rho ^{𝐫_1})[\rho (𝐫_2)\rho ^{𝐫_1}],$$
(7)
which can be rewritten exactly using Eq. (2) as

$$\rho (𝐫_1)=\rho ^{𝐫_1}+\rho ^{𝐫_1}𝑑𝐫_2c(r_{12};\rho ^{𝐫_1})[\rho (𝐫_2)\rho ^{𝐫_1}].$$
(8)
This is our final result, which we refer to as the HLR equation. We view this as an integral equation relating the hydrostatic density $`\rho ^{𝐫_1}`$ to the full density $`\rho (𝐫)`$ and solve it self-consistently for all $`𝐫_1`$. When $`\varphi (𝐫_1)`$ is known, we can immediately determine $`\rho ^{𝐫_1}`$ at each $`𝐫_1`$ from the local relation in Eq. (6), and then solve Eq. (8) by iteration for all $`𝐫_1`$ to determine the full density response $`\rho (𝐫)`$. Conversely, for a given equilibrium density distribution $`\rho (𝐫)`$ we can use HLR equation to determine the associated field $`\varphi (𝐫)`$. This inverse solution of Eq. (8) is particularly easy to carry out, since we can determine the local field at each $`𝐫_1`$ separately, without iteration. Accurate results have been obtained for many test cases with strong repulsive or attractive fields hlr ; katsov.weeks:incorporating .

This requires in particular expressions for $`\mu (\rho )`$ and for the direct correlation function $`c(r_{12};\rho )`$ of the uniform hard sphere fluid. In this paper we will use the Percus-Yevick (PY) py  approximation for $`c(r_{12};\rho )`$. This same function also arises from a self-consistent solution of the HLR equation, where the density response to a hard core field with $`R_v=\sigma `$ (equivalent to fixing a solvent particle at the origin) is related to the uniform fluid pair correlation function. Thus this self-consistent use of the HLR equation provides a physically suggestive way of deriving the PY result for $`c(r_{12};\rho )`$ hlr . The PY $`c(r_{12};\rho )`$ has a very simple analytical form and proves sufficiently accurate for our purposes here. Even better results can be found if one uses the very accurate expressions for the bulk $`c(r_{12};\rho )`$ and $`\mu (\rho )`$ as given by the GMSA theory gmsa  as inputs to the HLR equation.

### II.2 Relation to the PY approximation for a hard core solute

A spherical cavity acts like a hard core external field $`\varphi `$ that excludes the centers of all solvent molecules from the cavity region. We take the center of the cavity as the origin of our coordinate system, so that all distances are measured relative to the cavity center. Note that both the hydrostatic density $`\rho ^{𝐫_1}`$ from Eq. (5) and the full density response $`\rho (𝐫_1)`$ from Eq. (8) vanish whenever $`𝐫_1`$ is located in the cavity. This exact “hard core condition” comes out naturally from the theory, and does not have to be imposed by hand as in the GFM or the GMSA approaches.

To make contact with the PY approximation, recall that the cavity-solvent direct correlation function $`C(𝐫_1;\rho ^B,R_v)`$ for this system exactly satisfies

$$C(𝐫_1;\rho ^B,R_v)=𝑑𝐫_2\chi ^1(r_{12};\rho ^B)[\rho (𝐫_2)\rho ^B].$$
(9)
Thus $`C(𝐫_1)`$ is the function that replaces $`\beta \delta \varphi (𝐫_1)`$ so that the linear response equation (1) gives exact results when the full density change relative to the bulk is used on the right hand side. When this is compared to the HLR equation (7) for $`𝐫_1`$ outside the cavity region (where $`\rho ^{𝐫_1}=`$ $`\rho ^B`$ and $`\varphi =0`$) we see that the HLR equation predicts that $`C(𝐫_1)`$ vanishes. Thus for the HLR equation $`\rho (𝐫_1)`$ vanishes inside the cavity region and $`C(𝐫_1)`$ vanishes outside. This is the same as the PY approximation for the hard core cavity-solvent system hlr ; jdwannrev .

If $`R_v`$ is greater than $`\sigma /2`$, with $`\sigma `$ the solvent hard core diameter, then an equivalent exclusion is achieved by replacing the hard core external field by a hard core solute particle with (additive) diameter

$$\sigma _v2R_v\sigma .$$
(10)
From this it follows that if the PY approximation for the bulk $`c`$ or $`\chi ^1`$ is used, the density $`\rho (𝐫)`$ predicted by the HLR equation is identical to that given by the PY equation for the solute-solvent pair correlation function for a binary hard sphere mixture in the limit that the concentration of the solute species goes to zero HAB . Since an exact analytical solution of the PY equation for a binary HS mixture at arbitrary concentrations is known py-mixture-solution , we can take advantage of these results when computing the excess grand free energy.

This equality of solutions of the HLR equation and the PY mixture equation holds only for hard core cavity fields with radius $`R_v`$ $`\sigma /2`$ or $`\sigma _v0.`$ As discussed below in Sec. V, for tiny cavities with $`R_v`$ $`\sigma /2`$ the HLR equation can be solved directly and gives exact results for the density response if exact bulk correlation functions are used, and very accurate results when the PY approximation for the bulk $`c`$ is used. However, the corresponding PY mixture solutions in this range of $`R_v`$ (arrived at formally by taking $`\sigma _v`$ in Eq. (10) to be negative) are much less accurate. This inaccuracy arises from using the PY mixture solutions for negative $`\sigma _v`$. The direct solution of the PY cavity-solute equation for a tiny cavity, where a given approximation for the bulk $`c`$ is used along with the PY approximation that $`C(𝐫_1)`$ vanishes outside the cavity and $`\rho (𝐫_1)`$ vanishes inside, gives the same accurate results as the HLR equation. However, for more general external fields, the HLR equation is quite distinct from the PY approximation, and is generally more accurate; it has given good results for a wide range of fields hlr ; katsov.weeks:incorporating ; jdwannrev . This additional flexibility of the HLR equation will be required later in this paper when we discuss alternate density routes to the free energy.

## III Thermodynamic pathways to the free energy

In this section we discuss three different exact thermodynamic pathways for obtaining the excess free energy of inserting a cavity into a hard sphere fluid. The first two are well known, and the third describes a new density route that may have some computational advantages in other applications. We use the HLR equation to provide the needed structural information in all cases. We believe our calculation here represents the first use of a density route to obtain the excess free energy for this system. We then describe a simple but less accurate route to the free energy based on use of the partition function for the GFM.

### III.1 Compressibility route

In this route the excess free energy is determined by varying the chemical potential of the system while the external field $`\varphi (𝐫)`$ producing the cavity with radius $`R_v`$ remains constant. In the grand canonical ensemble the average number of particles $`N`$ is given by

$$\frac{\mathrm{\Omega }}{\mu }=N,$$
(11)
where $`\mathrm{\Omega }(\mu ,[\varphi ])`$ is the grand free energy. We can then calculate the free energy difference between the final state of interest and the trivial ideal gas state of zero density with $`\mu =\mathrm{}`$ and $`\mathrm{\Omega }=0`$ by integration:

$$\mathrm{\Omega }(\mu ^B,[\varphi ])=_{\mathrm{}}^{\mu ^B}𝑑\mu N=_{\mathrm{}}^{\mu ^B}𝑑\mu 𝑑𝐫\rho (𝐫;\mu ,[\varphi ]).$$
(12)
Then $`\mathrm{\Delta }\mathrm{\Omega }_v\mathrm{\Omega }(\mu ^B,[\varphi ])\mathrm{\Omega }(\mu ^B,[0]),`$ the desired excess grand free energy of the nonuniform fluid relative to the uniform bulk state, is given by

$$\mathrm{\Delta }\mathrm{\Omega }_v=_{\mathrm{}}^{\mu ^B}𝑑\mu 𝑑𝐫_1\{\rho (𝐫_1;\mu ,[\varphi ])\rho (\mu )\}.$$
(13)
As before, $`\rho (\mu )`$ gives the density of the uniform hard sphere solvent fluid as a function of the chemical potential.

Since $`\rho (𝐫;\mu ,[\varphi ])`$ vanishes inside the cavity, Eq. (13) shows there is a term in the excess free energy proportional to the cavity volume $`v`$ given by

$$v_{\mathrm{}}^{\mu ^B}𝑑\mu \rho (\mu )=v_0^{\rho ^B}𝑑\rho \frac{d\mu }{d\rho }\rho =vp^B,$$
(14)
on using the thermodynamic relation $`\rho (\mu /\rho )_T=(p/\rho )_T.`$ This exact leading order term for large $`v`$ is determined when using the compressibility route so that $`p^B`$ is the uniform fluid pressure calculated by the compressibility route hendersonSPT .

The term in curly brackets in Eq. (13) can be rewritten in a more convenient form for calculations by using the inverse relation to Eq. (9) for a general chemical potential $`\mu `$:

$$\rho (𝐫_1;\mu ,[\varphi ])\rho (\mu )=𝑑𝐫_2\chi (r_{12};\mu )C(𝐫_2;\rho (\mu ),R_v),$$
(15)
where $`\chi (r_{12};\mu )`$ $`\rho \delta (𝐫_1𝐫_2)+\rho ^2[g(r_{12})1]`$ is the usual linear response function of the uniform solvent fluid and $`g(r)`$ is the radial distribution function. Substituting into Eq. (13) and carrying out the integration over $`𝐫_1,`$ we have the formally exact result pdt

$`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ $`=`$ $`\beta {\displaystyle _{\mathrm{}}^{\mu ^B}}𝑑\mu \widehat{\chi }(0,\mu ){\displaystyle 𝑑𝐫_2C(𝐫_2;\rho (\mu ),R_v)}`$ (16)
$`=`$ $`{\displaystyle _0^{\rho ^B}}𝑑\rho \widehat{C}(0;\rho ,R_v).`$
Here $`\widehat{\chi }(0,\mu )`$ is the $`k=0`$ value of the Fourier transform of $`\chi `$, with a similar definition for $`\widehat{C}(0;\rho ,R_v)`$. In the last equality we used the uniform fluid compressibility relation $`\beta \widehat{\chi }(0,\mu )=d\rho (\mu )/d\mu `$ to change variables to an integration over density. We will explicitly solve the HLR equation for $`R_v\sigma /2`$ in Sec. V below, and from the equivalence between the HLR equation and PY mixture equation for $`R_v\sigma /2`$, we can use the exact solution of the PY mixture equation to obtain $`\widehat{C}`$ at larger $`R_v`$. Thus we can analytically carry out the integration in Eq. (16) for all $`R_v`$.

### III.2 Virial route

We now consider a different thermodynamic pathway, which was first used in scaled particle theory spt0 ; spt . Here we keep the chemical potential fixed at $`\mu ^B`$ and vary the range of the external hard core field by a scaling parameter $`\lambda `$, defining $`\varphi _\lambda (𝐫)\varphi (𝐫/\lambda )`$. For the hard core cavity field of interest here, as $`\lambda `$ is varied from $`0`$ to $`1`$ the radius of the exclusion zone then varies from $`0`$ to $`R_v`$. Since the density is generally related to the external field in the grand ensemble by

$$\frac{\delta \mathrm{\Omega }}{\delta \varphi (𝐫)}=\rho (𝐫),$$
(17)
the desired free energy difference is given by integration:

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=_0^1𝑑\lambda 𝑑𝐫\rho _\lambda (𝐫)\frac{\beta \varphi _\lambda (𝐫)}{\lambda }.$$
(18)
Here $`\rho _\lambda (𝐫)\rho (𝐫;\mu ^B,[\varphi _\lambda ]).`$ This formula is quite general and holds for any $`\lambda `$-dependent potential that vanishes for $`\lambda =0.`$ By exploiting special properties of the scaled hard core potential (the derivative of the Boltzmann factor of a hard core potential is a delta function) it is easy to show that Eq. (18) can be exactly rewritten as

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=4\pi R_v^3_0^1𝑑\lambda \lambda ^2\rho _\lambda (\lambda R_v).$$
(19)
Here $`\rho _\lambda (\lambda R_v)`$ is the contact density at the surface of the scaled exclusion zone with radius $`\lambda R_v`$. As in the compressibility route, we can analytically carry out the integration in the virial route to obtain solvation free energies for cavities for all $`R_v`$. The equivalent PY solution for binary hard sphere mixtures is used for the contact densities for all $`\lambda R_v`$’s larger than $`\sigma /2`$, while the explicit solution of the HLR equation is used for the $`\lambda R_v`$’s smaller than $`\sigma /2`$.

### III.3 Density routes

In addition to these particular pathways, we can also imagine directly changing the equilibrium density from $`\rho ^B`$ to $`\rho (𝐫)`$ over some convenient pathway specified by a coupling parameter $`\lambda ,`$ while taking account of the associated changes in $`\mathrm{\Omega }`$ and $`\varphi (𝐫)`$. Integrating Eq. (18) by parts to make $`\rho _\lambda (𝐫)`$ explicitly the controlling variable, we have exactly

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=𝑑𝐫\rho (𝐫)\varphi (𝐫)_0^1𝑑\lambda 𝑑𝐫\varphi _\lambda (𝐫)\frac{\rho _\lambda (𝐫)}{\lambda }.$$
(20)
Here $`\varphi _\lambda (𝐫)`$ is the external field consistent with the specified density profile, so that $`\rho (𝐫;\mu ^B,[\varphi _\lambda ])=\rho _\lambda (𝐫).`$ For a given density field $`\rho _\lambda (𝐫)`$, the HLR equation (8) can be solved inversely to obtain the associated hydrostatic density field $`\rho _\lambda ^𝐫`$. Using Eq. (6), $`\varphi _\lambda (𝐫)`$ at each $`𝐫`$ is locally related to $`\rho _\lambda ^𝐫`$ through $`\mu (\rho ).`$ Here we used the accurate Carnahan-Starling expression carnahan.starling:equation  for $`\mu (\rho )`$.

Most workers have considered a simple linear density path where

$$\rho _\lambda (𝐫)=\rho ^B+\lambda [\rho (𝐫)\rho ^B].$$
(21)
This has some theoretical advantages since $`\rho _\lambda (𝐫)/\lambda `$ is independent of $`\lambda `$ and has been successfully used in numerical calculations of the surface tension of the liquid-vapor interface katsov.weeks:incorporating . However, when $`\varphi (𝐫)`$ has a hard core (or is strongly repulsive), then the region near $`\lambda =1`$ in the $`\lambda `$-integration in Eq. (20) must be treated carefully, since for $`𝐫`$ in the hard core region $`\rho _\lambda (𝐫)/\lambda `$ is constant, while $`\varphi _\lambda (𝐫)`$ must tend to infinity as $`\lambda 1.`$ Although the singularity in the potential is integrable (for a hard core potential the divergent term in $`\beta \varphi _\lambda `$ goes as $`\mathrm{ln}(1\lambda )`$ and could be treated separately), in any case large contributions to the integral arise from a relatively small interval near $`\lambda =1`$. This could cause problems in a numerical integration.

To illustrate the computational advantages and flexibility that different pathways can provide, we introduce here a new density route that removes this possible difficulty. We consider a path that is linear in the square root of the density:

$$\rho _\lambda ^{1/2}(𝐫)(\rho ^B)^{1/2}+\lambda [\rho ^{1/2}(𝐫)(\rho ^B)^{1/2}],$$
(22)
where $`\rho _\lambda ^{1/2}(𝐫)=\sqrt{\rho _\lambda (𝐫)}`$ , etc. For this pathway we have

$$\frac{\rho _\lambda (𝐫)}{\lambda }=2\rho _\lambda ^{1/2}(𝐫)\frac{\rho _\lambda ^{1/2}(𝐫)}{\lambda }.$$
(23)
Both factors on the right side of Eq. (23) are easy to determine from Eq. (22). The numerical integration in Eq. (20) can now be carried out straightforwardly since the $`\rho _\lambda ^{1/2}(𝐫)`$ factor in Eq. (23) will cause $`\rho _\lambda (𝐫)/\lambda `$ to tend to zero exponentially fast wherever $`\varphi _\lambda (𝐫)`$ becomes large. Results using this path are reported below. Other paths implementing this idea exist and we have not tried to make an optimal choice.

### III.4 Gaussian Field Model

Finally we consider an alternative approach, the Gaussian field model (GFM) gft , that for hard core fields has many common elements with the HLR method. The GFM describes density fluctuations in a uniform fluid with average density $`\rho `$ by an effective quadratic Hamiltonian

$$^B=\frac{k_BT}{2}𝑑𝐫_1𝑑𝐫_2\delta \stackrel{ˇ}{\rho }(𝐫_1)\chi ^1(𝐫_{12};\rho )\delta \stackrel{ˇ}{\rho }(𝐫_2),$$
(24)
where $`\delta \stackrel{ˇ}{\rho }(𝐫)=`$ $`\stackrel{ˇ}{\rho }(𝐫)\rho `$ with $`\stackrel{ˇ}{\rho }(𝐫)`$ the microscopic density. The partition function for a system in an external field $`\varphi (𝐫)=\varphi _0(𝐫)+\varphi _1(𝐫),`$ with $`\varphi _0`$ a hard core field producing a cavity of radius $`R_v`$ and $`\varphi _1`$ a weaker perturbation, is then assumed to be given by gft

$`\mathrm{\Xi }_v`$ $`=`$ $`{\displaystyle 𝒟\stackrel{ˇ}{\rho }(𝐫)\left\{\underset{𝐫v}{}\delta [\stackrel{ˇ}{\rho }(𝐫)]\right\}}`$ (25)
$`\times \mathrm{exp}[\beta ^B+\beta {\displaystyle 𝑑r\stackrel{ˇ}{\rho }(𝐫)\varphi _1(𝐫)}].`$
The product of delta functions imposes the constraint that the density vanish inside the cavity. Inserting a Fourier representation for the $`\delta `$-functions and formally integrating $`\stackrel{ˇ}{\rho }(𝐫)`$ from $`\mathrm{}`$ to $`\mathrm{}`$ yields a Gaussian approximation for the partition function, as discussed below.

Moreover, using the same approximations, by functionally differentiating $`\mathrm{\Xi }_v`$ with respect to the field, one obtains the nonuniform singlet density in the GFM. In the case of a pure hard core field with $`\varphi _1=0,`$ the density response to a cavity with radius $`R_v`$ is given by

$$\rho (𝐫_1)=\rho ^B\rho ^B_v𝑑𝐫_2_v𝑑𝐫_3\chi (r_{12};\mu ^B)\chi _{in}^1(𝐫_2,𝐫_3).$$
(26)
The integrations are restricted to the cavity region, as indicated by the subscript $`v`$ on the integral symbols. Here $`\chi _{in}^1`$ is the inverse of the restricted linear response function $`\chi _{in}(r_{12};\mu ^B)`$, which equals $`\chi (r_{12};\mu ^B)`$ if both $`𝐫_1`$ and $`𝐫_2`$ are in the cavity region and equals zero otherwise. Thus $`\chi _{in}^1`$ is nonzero only inside the cavity and satisfies

$$_v𝑑𝐫_2\chi (r_{12};\mu ^B)\chi _{in}^1(𝐫_2,𝐫_3)=\delta (𝐫_1𝐫_3),$$
(27)
when both $`𝐫_1`$ and $`𝐫_3`$ are in the cavity region. Comparing Eq. (26) to Eq. (15), one can identify the cavity-solvent direct correlation function in the GFM as

$$C(𝐫_2)=\rho ^B_v𝑑𝐫_3\chi _{in}^1(𝐫_2,𝐫_3).$$
(28)
By properties of $`\chi _{in}^1,`$ the GFM $`C(𝐫)`$ vanishes outside the cavity region and $`\rho (𝐫)`$ in Eq. (26) vanishes inside. Thus the GFM gives exactly the same solution for the density response to a hard core external field as the PY or the HLR equations. (In the more general case where there is an additional perturbation potential $`\varphi _1`$, the various approaches differ. The GFM can be shown to treat the softer tail using the mean spherical approximation, which is different from and generally less accurate than the hydrostatic shift used in the HLR equation.)

Thus for cavities or hard core solutes all structurally based routes to the excess free energy will give the same results when using the GFM or the HLR equation. In addition, the GFM partition function also provides a direct and very simple route to the free energy gft . However this route is inherently approximate because Eq. (25) is not really a free energy functional for the whole configuration space, but rather a restricted one describing only the space outside of the specified cavity region. This functional may legitimately describe subsequent small perturbations of $`\varphi _1`$ outside of the cavity, but it does not contain enough information about the functional dependence on the cavity volume in the first place. Moreover the approximations made in evaluating the GFM partition function do not build in the fact that in grand canonical ensemble, the thermodynamic properties should depend on $`\mu \varphi `$ rather than on $`\mu `$ and $`\varphi `$ individually. Thus it is also not consistent with the free energy prediction from the compressibility route, which integrates over states at different chemical potentials but with a fixed hard core always present.

Evaluating the Gaussian integrals in Eq. (25), the excess grand free energy arising from a cavity with radius $`R_v`$ is given by

$`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ $`=`$ $`\mathrm{log}\mathrm{\Xi }_v/\mathrm{\Xi }^B`$ (29)
$`=`$ $`{\displaystyle \frac{1}{2}}\rho ^B\widehat{C}(0;\rho ^B,R_v)+\mathrm{log}(det\chi _{in}).`$
Here $`\mathrm{\Xi }_v`$ denotes the partition function with no particles in the cavity region and $`\mathrm{\Xi }^B`$ is the uniform bulk partition function. The term involving $`\widehat{C}`$ in Eq. (29) arises from the Gaussian integration of $`\stackrel{ˇ}{\rho }`$ in Eq. (25) and Eq. (28). When compared with the exact Eq. (16) from the compressibility route, we see the GFM effectively approximates $`\widehat{C}(0;\rho ,R_v)`$ for intermediate density values by $`\left(\rho /\rho ^B\right)\widehat{C}(0;\rho ^B,R_v)`$. This free energy contribution has a form similar to a harmonic oscillator with $`\left(\rho /\rho ^B\right)\widehat{C}(0;\rho ^B,R_v)`$ analogous to the restoring force. The second term is a result of the reduction of the configuration space.

An alternate perspective considers the average probability $`P_v(N)`$ of finding $`N`$ particles in the volume $`v`$ with radius $`R_v`$. The probability of inserting a cavity is thus $`P_v(0)`$. Then a formally exact expression for the excess free energy $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ is

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=\mathrm{log}P_v(0)=\mathrm{log}\frac{\mathrm{\Xi }_v[0]}{_{N=0}^{N_{max}}\mathrm{\Xi }_v[N]}.$$
(30)
Here $`\mathrm{\Xi }_v[N]`$ is the constrained partition function when $`N`$ particles are in the specified volume. This formula has been successfully used in the information theory approach developed by Hummer, Pratt and coworkers hummer.garde.ea:information .

If the GFM is used to approximate the partition functions in Eq. (30) by replacing the product of $`\delta `$-functions in Eq. (25) by the single average constraint $`\delta [_v𝑑𝐫\widehat{\rho }(𝐫)N]`$, one arrives at a Gaussian approximation hummer.garde.ea:information ; lcw  for $`\mathrm{\Xi }_v[N].`$ This “discrete” approximation for $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ based on this use of the GFM is

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=\mathrm{log}\frac{e^{\overline{N}^2/2\chi _v}}{_{N=0}^{N_{max}}e^{(N\overline{N})^2/2\chi _v}}$$
(31)
where

$$\chi _v_v𝑑𝐫_1_v𝑑𝐫_2\chi (|𝐫_1𝐫_2|),$$
(32)
and $`\overline{N}=\rho ^Bv`$.

We determined the solvation free energy $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ for the GFM using both the continuum version, Eq. (29), and the discrete version, Eq. (31). For the uniform $`\chi `$ we used the PY result. To estimate $`\mathrm{log}det\chi _{in}`$ we expanded $`\chi _{in}`$ in the volume $`v`$ using two orthogonal basis functions. A single constant basis function was used in Ref. lcw . This is exact for tiny cavities with $`R_v\sigma /2`$, as can be seen using Eq. (44) below. We chose one basis function to be constant. The other was taken to be $`j_0(R_vr/\pi )`$, the zeroth order spherical Bessel function with its first node fixed at $`r=R_v`$, but made orthogonal to the first (constant) basis function. The second basis function is thus a linear combination of $`j_0(R_vr/\pi )`$ and a constant. This was introduced to test the accuracy of the one basis function approximation previously used and hopefully will give improved results for larger $`R_v`$.

## IV Results for larger cavities with $`R_v>\sigma /2`$

We now discuss the solvation free energies given by the various pathways for a cavity with $`R_v>\sigma /2`$, equivalent to a physically realizable hard core solute particle with diameter $`\sigma _v>0.`$ (Results for tiny cavities with $`R_v\sigma /2`$ are discussed in Sec. V below.) We use the simplest version of the theory, where the PY approximation is used for the uniform fluid correlation functions. Fig. (1) give the solvation free energy $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ from the different pathways as a function of the packing fraction $`\eta =\pi \rho ^B\sigma ^3/6`$ for $`R_v/\sigma =1,\mathrm{\hspace{0.17em}1.5},\mathrm{\hspace{0.17em}1.75},`$ where the results can be compared to computer simulations gaus-sim  of Crooks and Chandler. Note that the volume of a spherical cavity with $`R_v=1.75\sigma `$ is over 42 times greater than that of a solvent particle. For $`R_v=\sigma `$ the results also give the excess chemical potential as a function of density for the uniform hard sphere fluid.

As discussed above, we can obtain analytical expressions for $`\mathrm{\Delta }\mathrm{\Omega }_v`$ for both the compressibility and virial routes. The compressibility route gives

$`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ $`=`$ $`{\displaystyle \frac{\eta (2+7\eta 11\eta ^2)}{2(1\eta )^3}}\mathrm{log}(1\eta )+{\displaystyle \frac{18\eta ^3}{(1\eta )^3}}{\displaystyle \frac{R_v}{\sigma }}`$ (33)
$`{\displaystyle \frac{18\eta ^2(1+\eta )}{(1\eta )^3}}{\displaystyle \frac{R_{v}^{}{}_{}{}^{2}}{\sigma ^2}}+{\displaystyle \frac{8\eta (1+\eta +\eta ^2)}{(1\eta )^3}}{\displaystyle \frac{R_{v}^{}{}_{}{}^{3}}{\sigma ^3}},`$
and the virial route gives

$`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ $`=`$ $`{\displaystyle \frac{\eta (2+7\eta 5\eta ^2)}{2(1\eta )^3}}\mathrm{log}(1\eta )`$ (34)
$`{\displaystyle \frac{18\eta ^2(1\eta )}{(1\eta )^3}}{\displaystyle \frac{R_{v}^{}{}_{}{}^{2}}{\sigma ^2}}+{\displaystyle \frac{8\eta (1+\eta 2\eta ^2)}{(1\eta )^3}}{\displaystyle \frac{R_{v}^{}{}_{}{}^{3}}{\sigma ^3}}.`$
Results for the density route and for the GFM are computed numerically.

In this range of $`R_v`$ there is good agreement except at the highest densities between the compressibility, virial and density routes, with best results overall arising from the compressibility route. The direct GFM predictions in Fig. (2) from the partition function are less satisfactory. Both the discrete and the continuum versions of the GFM give results that approach zero incorrectly as $`\rho ^B0`$, and the continuum values are consistently too large at high density while the discrete values are too small. The discrete version of the GFM uses a Gaussian approximation for the constrained partition functions and gives less accurate results than could be obtained from a fit to accurate values of $`N`$ and $`N^2`$ as in the information theory approach hummer.garde.ea:information .

As $`R_v\mathrm{}`$, the surface of the cavity approaches that of a planar wall. As shown in Eq. (14) there is a diverging term in the excess free energy given by the cavity volume $`v=4\pi R_v^3/3`$ times the bulk pressure $`p^B`$, and the more interesting quantity to calculate is the surface term $`\gamma _v`$, given by

$$\beta \gamma _v=\frac{\beta \mathrm{\Delta }\mathrm{\Omega }_v\beta p^Bv}{4\pi R_v^2}.$$
(35)
The surface tension of the planar wall is then $`\gamma _{\mathrm{}}`$. In the present case, both the compressibility and the virial routes give analytical expressions for $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ which depend on $`R_v`$ as a polynomial: $`a_0+a_1R_v/\sigma +a_2R_v^2/\sigma ^2+a_3R_v^3/\sigma ^3`$. The coefficient $`a_3`$ thus gives another route to the bulk pressure on taking the wall limit. The $`p^B`$ used in Eq. (35) has to agree with the prediction from the $`a_3`$ so that $`\gamma _{\mathrm{}}`$ is finite. As discussed earlier, the compressibility $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ yields the same bulk pressure as given by the accurate uniform fluid PY compressibility equation of state. The virial route does not automatically build in this consistency, and the pressure predicted from the coefficient $`a_3`$ is less accurate than the uniform fluid PY virial equation of state.

The coefficient of the quadratic term then gives the surface tension $`\beta \gamma _{\mathrm{}}=a_2/4\pi \sigma ^2`$. Using the compressibility route we find

$$4\pi \beta \gamma _{\mathrm{}}\sigma ^2=\frac{18\eta ^2(1+\eta )}{(1\eta )^3},$$
(36)
while the virial route gives

$$4\pi \beta \gamma _{\mathrm{}}\sigma ^2=\frac{18\eta ^2(1\eta )}{(1\eta )^3}.$$
(37)
The $`\gamma _{\mathrm{}}`$ obtained by the compressibility route coincides with that given by scaled particle theory spt0 ; spt ; j-henderson-simulation . We obtained $`\gamma _{\mathrm{}}`$ numerically for the density route.

These results can be compared to the quasi-exact formula d-henderson-formula

$$4\pi \beta \gamma _{\mathrm{}}\sigma ^2=\frac{18\eta ^2(1+\frac{44}{35}\eta \frac{4}{5}\eta ^2)}{(1\eta )^3},$$
(38)
which fits simulation data j-henderson-simulation  and imposes the known first and second surface virial coefficients bellemans . As shown in Fig. (3) the compressibility and density routes give excellent results, while the virial route is much less satisfactory.

This can be understood since the virial route uses only the contact densities at the fixed bulk density. The HLR equation is least accurate for the contact density at high bulk density and for large $`R_v`$, while the density response away from the solute is more accurate. On the other hand, the compressibility and density routes make use of the density response at all distances and over a range of densities from low density to the final $`\rho (𝐫)`$ where the HLR equation is more accurate.

## V A special regime: tiny cavities

### V.1 Exact results

The density response to a tiny cavity with $`R_v\sigma /2`$ is especially simple, since the center of at most one hard core solvent particle can lie anywhere within such a region spt0 ; spt . This fact allows one to determine exactly both the density response to a tiny cavity and the solvation free energy. Here we will compare these exact results to the predictions of the HLR and GFM approaches

As in Eq. (30), the solvation free energy is directly related to the average probability that no particle are in the cavity region

$$e^{\beta \mathrm{\Delta }\mathrm{\Omega }_v}=P_v(0),$$
(39)
where $`P_v(N)`$ is the probability of finding $`N`$ particles simultaneously in the region with volume $`v=4\pi R_v^3/3`$. When $`R_v\sigma /2`$ , the region can hold no more than one solvent particle, so that

$$P_v(0)+P_v(1)=1,$$
(40)
and

$$\rho ^Bv=N_v=P_v(1).$$
(41)
Substituting these into Eq. (39) we thus find the exact result spt0 ; spt

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=\mathrm{log}(1\rho ^Bv).$$
(42)
This argument can be extended to show that the exact density response to a tiny cavity is ReissC74 :

$$\rho (𝐫)=\frac{\rho ^B[1_v𝑑𝐫^{}\rho ^Bg(|𝐫^{}𝐫|)]}{(1\rho ^Bv)}$$
(43)
for $`𝐫`$ outside $`v`$, with $`\rho (𝐫)=0`$ for $`𝐫`$ inside. Here $`g(r)`$ is the exact radial distribution function for the uniform solvent fluid. Note that the contact density $`\rho (R_v)=\rho ^B/(1\rho ^Bv)`$ is exactly determined independent of the details of $`g(r)`$, since the corresponding $`g(|𝐫^{}𝐫|)`$ in Eq. (43) vanishes for all $`𝐫^{}`$ inside $`v`$. This result is valid as long as the inserted region $`v`$ can hold no more than one solvent particle, so Eq. (43) also holds for solvents with a hard core pair potential plus a softer tail.

### V.2 Structural predictions of the GFM and HLR methods

Now let us examine the GFM result in Eq. (26) for the case of a density response to a tiny cavity. Since $`g(|𝐫_1𝐫_2|)=0`$ when both $`𝐫_1`$ and $`𝐫_2`$ are in $`v`$, then $`\chi _{in}(𝐫_1,𝐫_2)=\rho ^B\delta (𝐫_1𝐫_2)(\rho ^B)^2`$ inside $`v`$. It is easy to see from Eq. (27) that the inverse function $`\chi _{in}^1`$ then has the simple form katsov:point

$$\chi _{in}^1(𝐫_1,𝐫_2)=\frac{\delta (𝐫_1𝐫_2)}{\rho ^B}+\frac{1}{1\rho ^Bv}.$$
(44)
When Eq. (44) is inserted into Eq. (26) to obtain the GFM density response to a tiny cavity, we recover the exact expression for $`\rho (𝐫)`$ given in Eq. (43), provided that the exact uniform fluid $`g(r)`$ or $`c(r)`$ is used. If approximate (say PY) results are used to describe the uniform fluid response functions then strictly speaking the GFM and HLR predictions for $`\rho (𝐫)`$ will not be exact for all $`𝐫.`$ However, the contact density $`\rho (R_v)`$ is exact in any case, since, as noted earlier, this requires only that the approximate $`g(r)`$ vanish inside the cavity region. Because of the equivalence between the structural predictions of the GFM and the HLR equation, these same conclusions hold for the HLR equation. In particular, the density response outside a tiny cavity is exactly described by linear response theory about the uniform bulk system.

### V.3 Solvation free energy from thermodynamic pathways

Moreover, these theories give the exact result of Eq. (42) for the solvation free energy $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$, independent of possible errors in $`g(r)`$, for all structurally based thermodynamic pathways that use only tiny hard core fields. In particular, the virial route in Eq. (19) gives exact results for $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ since this requires only the (exact) contact values $`\rho _\lambda (\lambda R_v).`$ The compressibility route as written in Eq. (16) requires $`C(𝐫_1)`$, which from Eq. (28) depends only on the exact $`\chi _{in}^1`$ in Eq. (44). Both results require only that the bulk $`g(r)`$ vanish for $`r<\sigma ,`$ and are unaffected by any errors at larger $`r`$. This is in accord with our general supposition that particular thermodynamic pathways can be relatively insensitive to errors in the structural theory.

However the choice of pathway is important. Thus the density routes do not give $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ exactly even in the tiny cavity regime. This is because as $`\rho _\lambda (𝐫)`$ is varied, the corresponding $`\varphi _\lambda (𝐫)`$ in general is not a pure hard core field and spreads outside the cavity region. Neither the GFM nor the HLR theories can treat these softer and longer-ranged fields exactly even if exact uniform fluid correlation functions are used.

### V.4 Solvation free energy from the GFM partition function

The $`\beta \mathrm{\Delta }\mathrm{\Omega }_v`$ obtained directly by taking the logarithm of the partition function in Eq. (29) can be expressed analytically for tiny cavities as

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=\frac{1}{2}\mathrm{log}(1\rho ^Bv)+\frac{1}{2}\frac{\rho ^Bv}{1\rho ^Bv}.$$
(45)
As $`v0`$, this goes as $`(\rho ^Bv)^2/4`$, while the exact result from Eq. (42) goes as $`\rho ^Bv`$. This quadratic term arises from the assumption of Gaussian fluctuations, which breaks down in this limit. The alternate discrete version from Eq. (31) reduces to

$$\beta \mathrm{\Delta }\mathrm{\Omega }_v=\mathrm{log}\frac{e^{\overline{N}/[2(1\overline{N})]}}{e^{\overline{N}/[2(1\overline{N})]}+e^{(1\overline{N})/[2\overline{N}]}},$$
(46)
where $`\overline{N}\rho ^Bv.`$ This has the peculiar behavior as $`v0`$ that all derivatives vanish, and so is significantly in error in this regime. See Fig. 4 for comparison with the exact answer.

## VI Conclusion

We have discussed several different thermodynamic routes that can be used to determine the solvation free energy for inserting both small and large cavities into a hard sphere fluid. Generally accurate results are found by using the HLR equation to relate the densities and associated fields over the intermediate states of the different pathways. We also considered the GFM and showed that it gives results equivalent to the HLR equation for the density response induced by a rigid cavity. However the GFM cannot describe the softer external potentials and intermediate densities needed for the density routes and for more general thermodynamic pathways. Direct use of the approximate partition function of the GFM to determine the solvation the free energy of a cavity also gives less accurate results.

Best results using the HLR equation for the solvation free energy of a cavity are found from the compressibility and density routes. This can be understood since most states along these routes require the density response at intermediate densities and distances away from the cavity where the HLR equation is most accurate. The HLR equation can also be used for more general solutes with different shapes or longer ranged attractive interactions and in applications where other pathways may be more useful. Combined with an appropriate pathway it represents a versatile and computationally efficient method for determining both the structure and the thermodynamics of nonuniform fluids.

We thank Kirill Katsov, Michael Fisher, Jim Henderson, and Lawrence Pratt for helpful comments and Gavin Crooks for sending us the simulation data that is used in the plots of finite size cavities. This work was supported by the National Science Foundation through Grant CHE-0111104.

*

## Appendix A

In the “wall” limit where $`R_v\mathrm{}`$, the surface tension given by the compressibility route is determined from

$$\gamma _{\mathrm{}}=_0^{\rho ^B}𝑑\rho \frac{\mu }{\rho }_0^{\mathrm{}}𝑑z[\rho _w(z)\rho ],$$
(47)
with $`\rho _w(z)\rho (z+R_v)`$. The virial route gives

$$\gamma _{\mathrm{}}=k_BT\rho ^B_0^{\mathrm{}}𝑑R[G(R)G(\mathrm{})],$$
(48)
Here $`\rho ^BG(R)`$ is the contact value of the density response to a cavity of size $`R`$. In the case considered in this paper, the $`\gamma _{\mathrm{}}`$ from both the compressibility and the virial routes can be determined analytically. In general, when no analytical expressions are available, one may need to carry out these integrations numerically.

There also exists the exact virial sum rule $`\rho ^BG(\mathrm{})=\beta p^B`$ for a planar wall immersed in the hard sphere fluid henderson:statistical . However, this focuses on the structure at contact for the planar wall, which is the region where the HLR equation is least accurate. Thus this thermodynamic pathway for the bulk pressure gives relatively poor results.