# Depletion potential in hard sphere fluids

## Abstract

A versatile new approach for calculating the depletion potential in a hard sphere mixture is presented. This is valid for any number of components and for arbitrary densities. We describe two different routes to the depletion potential for the limit in which the density of one component goes to zero. Both routes can be implemented within density functional theory and simulation. We illustrate the approach by calculating the depletion potential for a big hard sphere in a fluid of small spheres near a planar hard wall. The density functional results are in excellent agreement with simulations.

When the separation between two big colloidal particles suspended in a fluid of small colloidal particles or non-adsorbing polymers is less than the diameter of the small ones exclusion or depletion of the latter occurs leading to anisotropy of the local pressure which gives rise to an attractive depletion force between the big particles. Asakura and Oosawa  first described this depletion mechanism, suggesting that it would drive phase separation in colloid-polymer mixtures. Using excluded volume arguments they calculated the force for two hard spheres of radius $`R_b`$ in a fluid of small hard spheres of radius $`R_s`$ and showed that this is attractive for all separations less than $`2R_s`$ and is zero for larger separations. The hard sphere model captures the essence of the depletion phenomenon and can be mimicked experimentally by suitable choices of colloidal solutions . Depletion forces are of fundamental statistical mechanical interest since they are purely entropic in origin and much theoretical  and simulation  effort has been devoted to their determination. Many recent experiments, using a wide variety of techniques, have investigated depletion forces for colloidal mixtures or for colloid-polymer mixtures . Much of the impetus for these studies stems from a need to understand how depletion determines phase separation and flocculation phenomena and this has stimulated a growing interest in ascertaining quantitative details of the depletion force. Several of these studies are concerned with a big particle near a planar wall but Dinsmore et al.  have employed video microscopy to monitor the position of a single big colloid immersed in a solution of small colloids inside a vesicle – a system resembling hard spheres confined inside a hard cavity. The experimental conditions in Ref. correspond to a small sphere packing fraction $`\eta _s=0.3`$, which is sufficiently high that one expects the depletion force to be substantially different from that given by the Asakura-Oosawa result; the latter is valid only in the limit $`\eta _s0`$. Thus, in order to model this and other experimental situations one requires a theory of depletion which is reliable at high packings and which can tackle rather general “confining” geometries. The latter could represent a planar wall, a wedge or cavity or, indeed, another big particle . Here we present such a theory for the depletion potential based on a density functional treatment (DFT) of a fluid mixture. Our approach is versatile  and avoids the limitations of the virial expansion  and of the Derjaguin approximation  and avoids many of the limitations inherent in integral-equation theories  of the depletion force. It has a distinct advantage over a direct implementation of DFT in that it does not require a minimization of the free energy in the presence of the big particle. It also suggests an alternative simulation procedure for calculating depletion potentials. We demonstrate the accuracy of our approach by comparing our results for the depletion potential between a big hard sphere and a planar hard wall with those of our new simulation and with those of Ref..

We consider a multicomponent mixture in which species $`i`$ ($`i=1,2,\mathrm{}\nu `$) has a chemical potential $`\mu _i`$ and is subject to an external potential $`V_i(𝐫)`$ giving rise to the equilibrium number density profiles $`\{\rho _i(𝐫)\}`$. The quantity of interest $`W_t(𝐫_i)\mathrm{\Omega }_{ti}(𝐫_i;\{\mu _i\},T)\mathrm{\Omega }_{ti}(\mathrm{};\{\mu _i\},T)`$ is the difference in grand potential $`\mathrm{\Omega }_{ti}`$ between a test particle of species $`i`$ located at position $`𝐫_i`$ and one located at $`𝐫_i=\mathrm{}`$, i.e., deep in the bulk, far from the object exerting the external potentials. Using the potential distribution theorem in the grand ensemble one can easily show :

$$W_t(𝐫_i)=k_BT\mathrm{ln}\left(\frac{\rho _i(𝐫_i)}{\rho _i(\mathrm{})}\right)V_i(𝐫_i)+V_i(\mathrm{}),$$
(1)   
where $`T`$ is the temperature. Equation (1) is a general result; it is valid for arbitrary interparticle forces. It is important for subsequent application to note that $`\rho _i(𝐫)`$ is determined solely by the external potentials so it has the symmetry dictated by these potentials and not the broken symmetry brought about by inserting the test particle. In order to obtain the depletion potential $`W`$, which pertains to a single (big) particle, $`ib`$, in the presence of a fixed object exerting a potential $`V_b(𝐫)`$ of finite range, we take the limit $`\mu _b\mathrm{}`$ so that the density of the big particles vanishes, with the other chemical potentials $`\{\mu _{ib}\}`$ fixed. It follows that for $`𝐫`$ outside the range of $`V_b(𝐫)`$

$$\beta W(𝐫)=\underset{\mu _b\mathrm{}}{lim}\mathrm{ln}\left(\frac{\rho _b(𝐫)}{\rho _b(\mathrm{})}\right)$$
(2)   
with $`\beta =(k_BT)^1`$. Although both $`\rho _b(𝐫)`$ and $`\rho _b(\mathrm{})`$ vanish in this limit their ratio is non-zero. Moreover, the density profiles of the other species $`\{\rho _{ib}(𝐫)\}`$ reduce to those in the $`\nu 1`$ component fluid. Thus, for a binary mixture of big and small ($`s`$) particles, $`\rho _s(𝐫)`$ reduces to the profile of a pure fluid of small particles. Equation (2) can be re-stated in a more familiar form as $`p(𝐫)/p(\mathrm{})=\mathrm{exp}(\beta W(𝐫))`$, where $`p(𝐫)`$ is the probability density of finding the big particle at a position $`𝐫`$ from the fixed object.

There are two distinct ways of implementing this route to the depletion potential. First one can calculate the density profile $`\rho _b(𝐫)`$ and, hence, $`W_t(𝐫)`$ in the mixture for decreasing concentration of species $`b`$. For sufficiently negative values of $`\mu _b`$, $`W_t(𝐫)`$ should approach its limiting value $`W(𝐫)`$. Such a procedure is straightforward to implement in any approximate DFT, especially for objects with planar or radial symmetry (e.g., planar walls, spherical cavities or a fixed spherical particle). Second one can attempt to proceed immediately to the limit of a single big particle, thereby obtaining $`W(𝐫)`$ directly without evaluation of $`\rho _b(𝐫)`$. Equation (1) can be re-expressed as a difference in one-body direct correlation functions 

$$\beta W_t(𝐫_i)=c_i^{(1)}(𝐫_i;[\{\rho _i\}])c_i^{(1)}(\mathrm{};[\{\rho _i\}])$$
(3)   
and the depletion potential can be expressed formally as

$$\beta W(𝐫)=c_b^{(1)}(𝐫;[\{\rho _{ib}\},\rho _b=0])c_b^{(1)}(\mathrm{};[\{\rho _{ib}\},\rho _b=0]).$$
(4)   
In order to implement this result an explicit prescription for $`c_b^{(1)}`$ is required. Within the context of DFT $`c_b^{(1)}(𝐫_i;[\{\rho _i\}])=\beta \delta _{ex}[\{\rho _i\}]/\delta \rho _i(𝐫_i)`$, where $`_{ex}[\{\rho _i\}]`$ is the excess (over the ideal gas) free energy functional of the inhomogeneous mixture . For certain classes of (approximate) functionals explicit results can be obtained for $`c_i^{(1)}`$ which permit the limit $`\mu _b\mathrm{}`$ to be taken . An important example is Rosenfeld’s Fundamental Measure functional for hard sphere mixtures ,

$$_{ex}[\{\rho _i\}]=d^3r\mathrm{\Phi }(\{n_\alpha (𝐫)\}),$$
(5)   
where the Helmholtz free energy density $`\mathrm{\Phi }`$ is a function of the set of weighted densities

$$n_\alpha (𝐫)=\underset{i=1}{\overset{\nu }{}}d^3r^{}\rho _i(𝐫^{})\omega _i^\alpha (𝐫𝐫^{})$$
(6)   
and explicit expressions are available  for the weight functions $`\omega _i^\alpha `$ and for $`\mathrm{\Phi }`$. In this theory the free energy of the homogeneous mixture is identical to that from Percus-Yevick or scaled-particle theory. Within this approach Eq. (3) reduces to

$$\beta W_t(𝐫)=d^3r^{}\underset{\alpha }{}\mathrm{\Psi }^\alpha (𝐫^{})\omega _b^\alpha (𝐫^{}𝐫)$$
(7)   
with $`ib`$ and

$$\mathrm{\Psi }^\alpha (𝐫^{})\left(\frac{\beta \mathrm{\Phi }(\{n_\alpha \})}{n_\alpha }\right)_𝐫^{}\left(\frac{\beta \mathrm{\Phi }(\{n_\alpha \})}{n_\alpha }\right)_{\mathrm{}},$$
(8)   
i.e., the grand potential difference consists of a sum of convolutions of the functions $`\mathrm{\Psi }^\alpha `$ and the weight functions $`\omega _b^\alpha `$ of the big particle. The index $`\alpha `$ labels $`4`$ scalar plus $`2`$ vector weights . Since the derivatives in Eq. (8) are evaluated at equilibrium, the weighted densities must be obtained from Eq. (6), once the density profiles $`\rho _i(𝐫)`$ are obtained by minimizing the free energy functional. However, for the binary mixture the limit $`\rho _b(𝐫)0`$ implies that $`n_\alpha (𝐫)`$ involves only $`\rho _s(𝐫)`$ which is given by the solution of the Euler-Lagrange equation for the pure small-sphere fluid. The depletion potential $`W(𝐫)`$ is then given by Eq. (7) with $`n_\alpha (𝐫)`$ determined in this way. Both schemes have been used to calculate the depletion potential for a big hard sphere, immersed in a sea of small spheres, near a planar hard wall. Before describing our results we remark upon two limiting cases of the DFT treatment.

(i) The present functional  reduces to the exact low density limit $`\beta _{ex}[\{\rho _i\}]=\frac{1}{2}_{i,j}d^3rd^3r^{}\rho _i(𝐫)\rho _j(𝐫^{})f_{ij}(𝐫𝐫^{})`$, where $`f_{ij}`$ is the Mayer bond between species $`i`$ and $`j`$, when all the densities are small and the depletion potential simplifies to

$$\beta W(𝐫)=d^3r^{}(\rho _s(𝐫^{})\rho _s(\mathrm{}))f_{bs}(𝐫𝐫^{})$$
(9)   
where, in the same limit, $`\rho _s(𝐫)=\rho _s(\mathrm{})\mathrm{exp}(\beta V_s(𝐫))`$. For hard spheres $`f_{bs}(𝐫𝐫^{})=\mathrm{\Theta }((R_s+R_b)|𝐫𝐫^{}|)`$ where $`\mathrm{\Theta }`$ is the Heaviside function and one can show that for a hard potential $`V_s(𝐫)`$, $`\beta W(𝐫)`$ is $`\rho _s(\mathrm{})`$ times the overlap volume of the exclusion sphere of the big particle with the hard wall, i.e., one recovers the Asakura-Oosawa result for the depletion potential .

(ii) In the limit where the two species have equal radii, i.e., $`sR_s/R_b=1`$, it is straightforward to show that $`c_b^{(1)}`$ reduces to the direct correlation function of a pure ($`s`$) fluid and that the depletion potential is given by $`\beta W(𝐫)=\mathrm{ln}(\rho _s(𝐫)/\rho _s(\mathrm{}))`$, which is the exact result .

Since the Rosenfeld functional is known to yield very accurate results  for the density profile of a pure fluid near a hard wall and for the radial distribution function, obtained by fixing a particle at the origin and calculating the one-body density profile, we have good reasons to expect that the depletion potential calculated from DFT for $`s1`$ should be rather accurate for all (fluid) densities $`\rho _s(\mathrm{})`$. On the other hand for small size ratios the accuracy of the mixture DFT is not known; we test this by making comparison with simulations.

In Fig. 1 we show $`W_t(z)`$ for the binary hard sphere mixture, with $`s=0.2`$, at a planar hard wall as calculated using the Rosenfeld functional for a fixed packing fraction $`\eta _s4\pi R_s^3\rho _s(\mathrm{})/3=0.2`$ of the small spheres and 4 different values of $`\eta _b4\pi R_b^3\rho _b(\mathrm{})/3`$. $`z`$ is the distance between the surface of the big sphere and the wall. For $`\eta _b=10^4`$, $`W_t(z)`$ is indistinguishable from the depletion potential $`W(z)`$ calculated using the second, direct, method. Moreover, at this concentration $`\rho _s(z)`$ is indistinguishable from the the profile of the pure fluid. As both functions converge fairly rapidly to their limits (the results for $`\eta _b=10^2`$ differ from those for $`10^4`$ by at most 1.6%) we conclude that this is a practicable method of calculating depletion potentials via DFT. We note that for small values of $`\eta _b`$, $`W_t(z)`$ (and $`\rho _s(z)`$) exhibit decaying oscillations whose period is $`2R_s`$, whereas for $`\eta _b5\times 10^2`$ additional, longer period oscillations develop. The latter reflect the ordering of the big spheres which becomes pronounced at high packings. Figure 2 shows the depletion potential $`W(z)`$ calculated for a size ratio $`s=0.2`$ and $`\eta _s=0.3`$ alongside the result from Ref. . The latter was obtained from canonical Monte Carlo simulations of the depletion force, as given by an integral of the local density of small spheres in contact with the large sphere . Also plotted in Fig. 2 are our own simulation results, based on the formula $`p(z)\mathrm{exp}(\beta W(z))`$. As direct measurement of the probability density $`p(z)`$ of finding the big particle at a distance $`z`$ from the wall yields poor statistics, we measured the probability ratio $`p(z)/p(z+\mathrm{\Delta }z)`$ using umbrella sampling  in a grand canonical Monte Carlo simulation of the small spheres, i.e., we determined $`\beta (W(z+\mathrm{\Delta }z)W(z))`$, where $`\mathrm{\Delta }z`$ is a step length, and we set $`W(0)`$ equal to the DFT value. Our DFT results are in excellent agreement with those of both simulations: the height of the depletion barrier (the maximum value of $`W`$ minus the contact value $`W(0)`$) is given almost exactly and the subsequent extrema are closely reproduced by the theory . That two independent sets of simulation results, based on totally different routes to the depletion potential, agree so closely is pleasing and attests to the accuracy with which these potentials can be calculated. We find a similar high level of agreement between DFT and the results of Ref.  for $`\eta _s=0.1,0.2`$ and the smaller size ratio $`s=0.1`$.

The oscillatory behavior of $`W(z)`$ warrants further discussion. From the general theory of the asymptotic decay of density profiles in mixtures with short-ranged interaction potentials it is known that the profiles of all species near a wall (or a fixed particle) exerting a finite ranged external potential exhibit the same characteristic decay . Thus, for a binary mixture near a hard wall $`\rho _i(z)\rho _i(\mathrm{})A_{wi}\mathrm{exp}(a_0z)\mathrm{cos}(a_1z\mathrm{\Theta }_{wi})`$, $`z\mathrm{}`$ where $`i=b,s`$. The amplitude $`A_{wi}`$ and phase $`\mathrm{\Theta }_{wi}`$ depend on the particular species but $`a_0`$ and $`a_1`$ are common. In the limit $`\rho _b(\mathrm{})0`$ these are given by the solution of $`1\rho _s(\mathrm{})\widehat{c}_{ss}(a)=0`$, where $`\widehat{c}_{ss}(a)`$ is the Fourier transform of the two-body direct correlation function of the homogeneous pure fluid of density $`\rho _s(\mathrm{})`$. $`a_1`$ is the real and $`a_0`$ the imaginary part of the solution $`a`$ with the smallest imaginary part . We have confirmed from our numerical results that for $`\eta _b=10^4`$ and $`s=0.2`$, $`\rho _b(z)`$ and $`\rho _s(z)`$ exhibit the same value of $`a_0`$ and $`a_1`$. Moreover, for fixed $`\eta _s`$, both $`a_0`$ and $`a_1`$ are, as predicted, independent of the size ratio $`s`$. It follows from Eq. (2) that $`\beta W(z)A_{wb}\mathrm{exp}(a_0z)\mathrm{cos}(a_1z\mathrm{\Theta }_{wb})`$, $`z\mathrm{}`$, and apart from the amplitude and phase, the asymptotic decay of the depletion potential is also independent of $`s`$ . From a fit to the numerical data for $`\eta _s=0.2`$ we obtained $`a_02.1/(2R_s)`$ and $`a_15.2/(2R_s)`$. These values are close to those obtained from solving $`1\rho _s(\mathrm{})\widehat{c}_{ss}(a)=0`$ for the pure fluid where, consistent with the Rosenfeld DFT, $`\widehat{c}_{ss}`$ is given by Percus-Yevick theory. Moreover, and in keeping with earlier studies of bulk pair correlation functions , we find that the leading-order asymptotic result provides a good description of the intermediate as well as the longest-ranged behavior of the depletion potential. We note that the period of the oscillations $`2\pi /a_1`$ decreases whereas the decay length $`a_0^1`$ increases with increasing $`\eta _s`$.

The success of the comparisons between our DFT results and those of simulation lead us to expect that the Rosenfeld functional will provide an accurate description of depletion potentials for size ratios down to $`s=0.1`$ and for packing fractions $`\eta _s`$ up to $`0.3`$. (For smaller size ratios we expect the Percus-Yevick theory of hard sphere mixtures, upon which the functional is based, to become inaccurate and there is no reason to expect that any approximate mixture functional would reproduce the exact $`s0`$ results for the depletion potential which are obtained by making the Derjaguin approximation at the outset .) This implies that our procedure can be profitably employed for some of the more complex geometries investigated in experiments . We emphasize that our present procedure has a crucial advantage over the brute-force application of DFT to the calculation of depletion potentials. In the latter one would calculate $`W(𝐫)`$ either from the total free energy of the inhomogeneous fluid or from the local density of the small particles in contact with the big one . In contrast with the present scheme, both methods require much numerical effort to minimize the free energy functional since the original symmetry of the density profile $`\rho _s(𝐫)`$ is broken by the presence of the big particle. A likely limitation of our procedure is that one requires a reliable free energy functional for the (binary) mixture of big and small particles. Although accurate functionals are available for hard sphere mixtures  this is not the case for most other types of mixtures.

As a final remark we note that it is possible to take simulation data for $`\rho _s(𝐫)`$, computed in the absence of the big test particle, and insert these into the DFT expression Eq. (7). Although such a procedure is not self-consistent, in view of the high accuracy of the Rosenfeld functional it does provide an alternative means of calculating depletion potentials for complex geometries such as wedges where a direct simulation of the potential or force is extremely difficult.

We are grateful to R. Dickman for providing his simulation results for the depletion potential which we plot in Fig. 2. We benefited from conversations with R. Van Roij. This research was supported in part by the EPSRC under GR/L89013.