{'fraction_non_alphanumeric': 0.044679239839890375, 'fraction_numerical': 0.013865349248133858, 'mean_word_length': 4.370678803137407, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hard-sphere fluids in annular wedges: density distributions and depletion potentials.

We analyze the density distribution and the adsorption of solvent hard spheres in an annular slit formed by two large solute spheres or a large solute and a wall at close distances by means of fundamental measure density functional theory, anisotropic integral equations and simulations. We find that the main features of the density distribution in the slit are described by an effective, two--dimensional system of disks in the vicinity of a central obstacle. For large solute--solvent size ratios, the resulting depletion force has a straightforward geometrical interpretation which gives a precise"colloidal"limit for the depletion interaction. For intermediate size ratios 5...10 and high solvent packing fractions larger than 0.4, the explicit density functional results show a deep attractive well for the depletion potential at solute contact, possibly indicating demixing in a binary mixture at low solute and high solvent packing fraction.


I. INTRODUCTION
The equilibrium statistical theory of inhomogeneous fluids whose foundations were laid out in the 1960's [1,2] is a well-studied subject which has been developed since by a fruitful interplay between simulations and theory. For the simplest molecular model, the hard-body fluid mixture, this has led to the development of powerful, geometry-based density functionals [3,4,5,6,7,8] which accurately describe adsorption phenomena, phase transitions (such as freezing for pure hard spheres or demixing for entropic colloid-polymer mixtures) and molecular layering near obstacles. Such a class of density functionals has not been found yet for fluids with attractions, however, significant progess has been achieved in the description of bulk correlation functions and phase diagrams through the method of integral equations [9,10,11].
Strong inhomogeneities occur if fluids are confined on molecular scales, a topic which enjoys continuous interest [12,13]. For example, the packing of molecules between parallel walls leads to oscillating forces between them which can be measured experimentally and determined theoretically [14]. With the development of preparational techniques for colloids and their mixtures, the "molecular" length scale has been lifted to the range of several nm up to µm which even allows the direct observation of modulated density profiles through microscopy besides the measurement of resulting forces on the walls. Furthermore, in the colloidal domain the possibility to tailor the interparticle interactions to a certain degree allows the close realization of some of the favourite models for simple fluids fancied by theorists, such as e.g. hard spheres [15]. This has opened the route to quantitative comparisons between experiment, theory and simulations.
In an asymmetric colloidal mixture with small "solvent" particles and at least one species of larger solute particles, the phenomenon of solvent-mediated, effective interactions between the solute particles may give rise to various phase separation phenomena [16]. From a theoretical point of view, these effective interactions are interesting since they facilitate the description of mixtures in terms of an equivalent theory for one species interacting by an effective potential [17,18]. If the solvent particles possess hard (or at least steeply repulsive) cores, they are excluded from the region between two solute particles if the latter are separated by less than one solvent diameter. This gives rise to strong depletion forces whose understanding is crucial for the concept of an effective theory containing only solute degrees of freedom. Since the magnitude of the depletion force is directly linked to the solvent density distribution around the solutes, the quantitative investigation of the solvent confined between solute particles appears to be important. The confinement becomes rather extreme for large asymmetry between solute and solvent (see below).
In the following we want to concentrate on the idealized system of additive hard spheres and in particular on the effective interaction between two solute particles in solvent, i.e. the case of infinite dilution of solute particles in the colloidal mixture. The case of the interaction of one solute particle with a wall is a special case in that the radius of the other solute particle goes to infinity. For the case of hard spheres (solvent diameter σ = 2R s , solute radius R b such that the size ratio is α = R b /R s and R = R b + R s is the radius of the exclusion sphere for solvent centres around a solute sphere) it has been found that the theoretical description using "bulk" methods becomes increasingly inaccurate for size ratios α 5 and solvent densities ρ * s = ρ s σ 3 0.6, when compared to simulations. Here, the terminus "bulk" methods refers to methods which determine the bulk pair correlation function between the colloids, g bb (r), from which the depletion potential is obtained as βW = − ln g bb . Bulk inte-gral equation (IE) methods [19] and the "insertion" trick in density functional theory (DFT) [20,21,22,23] fall into this category. For larger size ratio, one would expect that the wellknown Derjaguin approximation [24,25,26] becomes accurate very quickly. Interestingly, however, both the mentioned bulk methods and simulations deviate significantly from the Derjaguin approximation (besides disagreeing with each other) for size ratios 10 and solvent densities ρ * s = 0.6 . . . 0.7 [26,27,28] when the surface-to-surface separation h between the solute particles is close to σ (h σ, i.e. at the onset of the depletion region where the solvent particles are "squeezed out" between the solute particles). One notices that for h ≤ σ the solute particles form an annular wedge with a sharp edge which restricts the solvent particles to quasi-2d motion (see Fig. 1). Taking this observation into account, the phenomenological analysis of Ref. [26] predicted that the solvent adsorption at this edge leads to a line contribution to the depletion potential which is proportional to the circumference of the circular egde and thus to R 1/2 . Such a term in the depletion potential causes a very slow approach to the Derjaguin limit for large solutes (the Derjaguin approximation describes the depletion potential essentially by volume and area terms of the overlap of the exclusion spheres pertaining to the solutes [25,26], see below). In this way, a generalized Derjaguin approximation serving as a new "colloidal" limit can be formulated which turns out to have a far more general meaning [29] than anticipated. Using the concept of morphological (morphometric) (morphometric) (morphometric) (morphometric) (morphometric) (morphometric) (morphometric) (morphometric) thermodynamics introduced in Ref. [30], one finds that the insertion free energy of two solute particles (and thus the depletion potential between them) only depends on the volume, surface area, and the integrated mean and Gaussian curvatures of the solvent accessible surface around the two solutes. The coefficients of these four terms depend on the solvent density but not on the specific type of surface. In this manner, the phenomenological line tension of Ref. [26] is related to the general coefficient of mean curvature of the hard-sphere fluid [29], and the Derjaguin approximation is equivalent to the morphometric analysis restricted to volume and surface area terms.
Since the depletion force in the hard-sphere system is directly linked to an integral over the solvent contact density on one solute (see Eqs. (1) and (2) below), one should be able to connect the morphometric approach to features of the solvent density profile in the annular wedge. This is a strong motivation for us to explicitly determine the wedge density profile.
We will do so by means of density functional theory and anisotropic integral equations, and compare it to results of Monte-Carlo (MC) simulations for selected parameters. For the intermediate densities ρ * s = 0.6 and 0.7 and various size ratios between solute and solvent results for the depletion force from such explicit DFT calculations were already presented in Ref. [29], and good agreement with the morphometric depletion force was obtained for size ratios between 5 and 40. Here, we will present a detailed analysis of the wedge density profiles and consider also higher densities. We will give a thorough comparison to recent results from MC simulations which have been obtained for the wall-solute interaction at a solvent density ρ * s = 0.764 (η s = (π/6)ρ s σ 3 = 0.4) and size ratios between 10 and 100 [31,32,33].
The paper is structured as follows. In Sec. II we introduce the basic notions of density functional and integral equation theory and present a short description of the Monte Carlo simulations employed here. In Sec. III A we analyze in detail the density profiles and the corresponding depletion forces for the wall-sphere geometry for the particular solvent density ρ * s = 0.764. This geometry permits explicit calculations up to solute-solvent size ratios 100 and a test of the "colloidal limit" of the morphometric depletion force. Sec. III B gives an analysis of the depletion force and potential in the sphere-sphere geometry (corresponding to the important case of the effective interaction in a dilute mixture) for intermediate size ratios 5 and 10 and higher solvent densities ρ * s = 0.8 . . . 0.9.

II. THEORY AND METHODS
We consider two scenarios (see Fig. 2): (a) One hard solute sphere immersed in the hard sphere solvent of density ρ s confined to a slit, created by two hard walls at distance L.
The surface-to-surface distance between one wall and the solute sphere is denoted by h, furthermore L ≫ h such that the correlations from one wall do not influence the solute interaction with the other wall. (b) Two hard solute spheres immersed in the bulk solvent spheres with surface-to-surface distance h. The solvent accessible surface is given by the dashed lines in Fig. 2, thereby one sees that for h ≤ σ the two solutes (or the solute and one wall) form an annular wedge in which the solvent adsorbs. The solvent density profile ρ(r) ≡ ρ(r || , z) depends only on z, the coordinate on the symmetry axis and the distance r || to the symmetry axis. In the latter case (b), when the first colloid is centered at the origin and the second one at z = 2R b + h, the total force f (h) on the first colloid is obtained by integrating the force −∇u bs (r) (u bs is the solute-solvent potential) between the solute and one solvent sphere over the density distribution ρ(r): Thus the total depletion force reduces to an integral over essentially the contact density on the exclusion sphere around the colloid. This follows from exp(−βu bs )β∇u bs = −rδ(|r| − R) and the observation that ρ(r) exp(βu bs ) is continuous across the exclusion sphere surface. In case (a), the force on the colloid is the negative of the excess force f w on the wall and the excess force is determined by the total force on the wall minus the contribution from the bulk solvent pressure p. Through the wall theorem, the latter is given by βp = ρ w where ρ w is the contact density of the solvent at a single wall. By an argument similar to the above one and putting the exclusion surface of the wall at z = 0, f w is determined as Thus, f w is equivalent to the excess adsorption at the wall which makes it somewhat easier to determine in MC simulations than f [32].

A. Density functional theory
The equilibrium solvent density profile ρ(r) ≡ ρ eq (r) can be determined directly from the basic equations of density functional theory. The grand potential functional is given by where F id and F ex denote the ideal and excess free energy functionals of the solvent. The solvent chemical potential is denoted by µ and the solute(s) and/or the walls define the external potential V ext . The ideal part of the free energy is given by with Λ denoting the de-Broglie wavelength. The equilibrium density profile ρ eq (r) for the solvent at chemical potential µ = β −1 ln(ρ s Λ 3 ) + µ ex (corresponding to the bulk density ρ s ) given by σ, and the radius of the exclusion sphere around a solute particle is given by is found by minimizing the grand potential in Eq. (3): For an explicit solution, it is necessary to specify the excess part of the free energy. Here we employ two functionals of fundamental measure type (FMT). These are the original Rosenfeld functional (FMT-RF) [3] and the White Bear functional (FMT-WB, with Mark I from Refs. [4,5] and Mark II from Ref. [6]), for another closely related variant see Ref. [7].
It has been demonstrated that FMT gives very precise density profiles also for high densities of the hard sphere fluid in various circumstances. Furthermore, the White Bear II functional possesses a high degree of self-consistency with regard to scaled-particle considerations [6].
However, we do not consider the tensor-weight modifications of these functionals which are necessary to obtain a correct description of the liquid-solid transition [34] and are of higher consistency in confining situations which reduce the dimensionality of the system ("dimensional crossover", see Ref. [35,36]). This might be an issue in some circumstances (see below).
Th explicit forms of the excess free energy are given in App. B. Numerically, due to the non-local nature of F ex , Eq. We remark that there have been earlier attempts to obtain explicit density distributions using DFT around two fixed solutes and thereby to extract depletion forces [37,38,39]. In Ref. [38] this was done using the simple Tarazona I functional for hard spheres [40], and considerations were limited to a small solvent packing fraction of 0.1 and solute-solvent size ratios of 5. In Ref. [39], a minimization of FMT-RF was carried out using a real-space technique for solvent densities up to ρ * s = 0.6 and size ratios 5. In that respect, the present technique is superior in that it allows us to present solutions for solvent densities up to ρ * s = 0.9 and size ratios up to 100.

B. Integral equations
The central objects within the theory of integral equations are the correlation functions on the one and two-particle level in the solvent. We are interested in the two-particle correlation functions in the presence of an external background potential (one fixed object, solute or wall), from which the explicit solvent density profiles around the two fixed objects (solute-solute or solute-wall) follow (see below). This is usually referred to as the method of inhomogeneous or anisotropic integral equations, first employed for Lennard-Jones fluids on solid substrates [41,42] and for hard-spheres (HS) fluids in contact with a single hard wall in Ref. [43]. Liquids confined between two parallel walls (a planar slit) have been extensively studied by Kjellander and Sarman [44,45,46].
For an arbitrary background potential V (r) which corresponds to an equilibrium background density profile for the solvent ρ V (r), the pair correlation function g ij (r, r 0 ) = h ij (r, r 0 ) + 1 describes the normalized probability to find a particle of species i at position r if another particle of species j is fixed at position r 0 . For our system, the species index is either s (solvent particle) or b (big solute particle). In the solute-solute case, V is given by the potential of one solute particle, whereas in the solute-wall case it is the potential of the wall. Then the equilibrium density profile discussed in the previous subsection is related to the pair correlation function through ρ(r) = ρ V (r) g sb (r, r 0 ) (r 0 specifies the position of the (other) solute particle). The depletion force follows then through Eqs. (1) and (2).
The corresponding direct correlation functions of second order c ij (r, r 0 ) are related to h ij through the inhomogeneous Ornstein-Zernike (OZ) equations In the dilute limit for the solute particles which we consider here, the background density ρ V,b for the solutes is zero, therefore the OZ equations reduce to The background density profile is linked to c ss through the Lovett-Mou-Buff-Wertheim equation [47]: A third set of equations is necessary to close the system of equations. The diagrammatic analysis of Ref. [1] provides the general form of this closure which reads: where u ij are the pair potentials in the solute-solvent mixture and b ij denote the bridge functions specified by a certain class of diagrams which, however, can not be resummed in a closed form. For practical applications, these bridge functions need to be specified in terms of h ij and c ij to arrive at a closed system of equations. Among the variety of empirical forms devised for this connection we mention those which, in our opinion, rely on somewhat more general arguments. These are the venerable hypernetted chain (HNC), Percus-Yevick (PY) closure and the mean spherical approximation (MSA) which can be derived by systematic diagrammatic arguments [48], and the reference HNC closure which introduces a suitable reference system for obtaining b ij with subsequent free energy minimization [49]. (For bulk properties of liquids, the self-consistent Ornstein-Zernike approximation (SCOZA) [9] and the hierarchical reference theory (HRT) [10] are very successful. They rely on a closure of MSA type but it appears to be difficult in generalizing them to inhomogeneous situations such as considered here.) In calculations, we considered the closures where γ ij = h ij − c ij . The Percus-Yevick (PY) approximation Eq. (11) is exactly solvable in the bulk case even for the muticomponent HS fluid [50]. However, beyond providing good qualitative behavior PY does not produce precise quantitative results in general, since it fails at contact, where the value of pair distribution function is too small. For bulk systems, PY generates a noticeable thermodynamic inconsistency, i.e. a discrepancy between different routes to the equation of state. To overcome these deficiencies a refined approximation to the closure, which interpolates between HNC (b ij = 0) and PY closures was suggested by Rogers and Young (RY) [51], see Eq. (12) where the ξ ij are adjustable parameters, which can be found from the thermodynamic consistency requirement [51]. An even more suitable for the asymmetric HS mixtures variant of the modified Verlet (MV) closure was suggested in Ref. [52], see Eq. (13). Here the parameters α ij are chosen to satisfy the exact relation between b ij (0) and the third virial coefficient at low densities. For the choice of the auxiliary parameters ξ ij and α ij in the present work, see App. C We have solved the set of equations (7)-(13) for the solute-wall case where we could employ similar numerical methods as in the numerical treatment of density functional theory.
(This is not possible for the solute-solute case, for an efficient method applicable in this case, see Ref. [53].) In the solute-wall case, the background density profile ρ V (r) ≡ ρ V (z) depends only on the distance z from the wall. The two-particle correlation functions depend on the z-coordinates of the two particles individually and the difference in the radial coordinates: . Thus, the Ornstein-Zernike equations (7) and (8)

C. Simulation Details
The Monte Carlo simulations for the wall-solute system were performed at fixed particle number and volume of the simulation box. The upper boundary of the box was given by a large hard sphere, the lower boundary by a planar hard wall (see hatched regions in Fig. 3).
All remaining boundaries were treated as periodic. We sampled the configuration space of the small spheres by standard single particle translational moves. In order to impose the asymptotic solvent bulk density ρ * s = 0.764 (η s = 0.4), the concentration of small spheres in the box was set such that the density at contact with the hard wall far away from the wedge (i. e. in the region marked by grey squares in Fig. 3, averaged over a depth 0.02 σ) settled toρ Wall = 4.88 within 1% error. This valueρ Wall = 4.88 was obtained from the density profile of hard spheres at a hard wall at the bulk packing fraction 0.4, calculated with FMT-WBII which is very accurate. System sizes ranged from 1800 particles for R b = 10 σ to 8000 particles for R b = 25 σ. In order to access the configurations inside the narrow part of the wedge with sufficient accuracy, large numbers of Monte Carlo sweeps were required.
We equilibrated the systems for 5 × 10 5 MC sweeps (i. e. attempted moves per particle) each. For data acquisition, we performed between 10 7 and 10 8 sweeps (depending on the parameters h and R b ) and averaged over 10 6 to 10 7 samples. Note that from the simulations we obtained only the wedge density profiles (see Figs. 5-7 below) but did not attempt to obtain the depletion force from Eq. (2) where one needs an excess adsorption integral at the wall. This would require considerably more sweeps [32,33], see also

A. Wall-solute interaction
In this section, we analyze case (a), the problem of a solute sphere close to a wall. Previous simulation work [31,32,33] has concentrated on the particular state point ρ * s = 0.764 (packing fraction η s = 0.4), therefore we present explicit results for the density profiles also for this state point to facilitate comparison.

Density profiles
In Fig. 4, we show DFT results for the solvent density profile between wall and solute for a solute-solvent size ratio of 20 and solvent-wall distance h = 0 (i.e. contact between wall and solute molecular surface; left panel) and h = σ (i.e. the end of the depletion region; right panel). There is strong adsorption at the apex of the annular wedge, given by the coordinates z = 0 and r || = r 0 = 2R(h − σ) − (h − σ) 2 , as reflected by the main peak.
Along the wall (z = 0), we observe strong structuring which is mainly dictated by packing considerations. The apex peak corresponds to a ring of solvent spheres centered at r || = r 0 (for h = 0) or just one sphere near r || = 0 (for h = σ). The second-highest peak in both show a suitable weighted average of the zeroth and first bin (∆z = 0.05), assuming linear interpolation. As seen in the right panel of Fig. 5, IE-MV gives a much lower density close to the apex of the wedge compared to the other methods. This is presumably due to the lower resolution in z-direction which is possible for IE (see App. C).
To further investigate the issue of the dominant packing mode in the annular wedge, we performed DFT calculations also for the additonal size ratios 10, 30, 50 and 100 for the same solvent density ρ * s = 0.764. Packing in r || -direction can be monitored conveniently by viewing the annular wedge close to the apex as a quasi-2d system (see Fig. 1 and the remarks in the Introduction). We define an effective, two-dimensional solvent density by where defines the surface of the exclusion sphere around the solute (and thus the radius-dependent width of the annular wedge). In Fig. 6 we show results for ρ 2d (r || ) for two configurations with slightly different solute-wall distances: h = 0.95 σ (left panel) and h = σ (right panel). In the latter case, one solvent sphere fits exactly between solute and wall at the apex, nevertheless ρ 2d vanishes there due to the vanishing slit width, z s (r || → 0) → 0. Note that this is also a consequence of an application of the potential distribution theorem [54] which states that the 3d density reaches a finite, maximum value of exp(βµ ex ) in a planar slit of vanishing width which is equivalent to a vanishing 2d density (µ ex is the excess chemical potential of the bulk solvent coupled to the slit). However, when the slit widens, the 2d density quickly rises [26]. For h = σ (right panel in Fig. 6), this leads to the picture of an effective 2d system with a small, soft and repulsive obstacle centered at r || = 0 which induces moderate layering in the 2d-"bulk" with an effective "bulk" density ρ 2d,s ≈ 0.7 . . . 0.8. For h = 0.95 σ (left panel in Fig. 6), the obstacle in the center has a radial extent of r 0 ≈ Rσ/10 which ranges from 0.75 σ (size ratio 10) to 2.25 σ (size ratio 100) and induces stronger layering in the effective 2d system. The oscillations occur around the same "bulk" density as in the case h = σ.
The simulation results show weaker oscillations than the DFT results. Since the first peak, e. g., occurs at wedge widths < 0.03 σ, we are testing the dimensional crossover properties of DFT-FMT from 3d to 2d with a sensitive probe. The difference between the results of the simulations and the DFT versions employed here are in line with previous investigations on the dimensional crossover properties of FMT [36]. There it has been found that the strict 2d limit of FMT-RF results in a somewhat peculiar (integrable) divergence in the hard disk direct correlation function c 2d (r) for r → 0 and overestimated peaks in the corresponding structure factor. The tensor weight modifications introduced in Ref. [34] result in better dimensional crossover properties [56] and might improve the DFT results close to the apex of the annular wedge. A more detailed investigation of tensor-weighted FMT with regard to the correlations in narrow planar slits and also in the annular wedges considered here is certainly of interest (although the algebra of App. B becomes much more extended in this case).
Our observations should be compared with the phenomenological modelling of Ref. [26] ("annular slit approximation"). There, the effective 2d system was approximated by an idealized system of hard disks, layering around a hard cavity of radius r 0 centered around r || = 0. The bulk density ρ 2d,s of this hard disk system was determined via the selfconsistency condition p 2d (ρ 2d,s ) = −2γ(ρ s ) where p 2d is the 2d pressure and γ is the surface tension of a hard wall immersed in a bulk solvent of density ρ s . Employing scaled particle theory [55], this yields ρ * 2d,s ≈ 0.66 for ρ * s = 0.764. This value is a bit smaller than the one inferred from the 2d density profiles of Fig. 6. The main difference between the idealized hard disk system of Ref. [26] and the effective 2d system as reflected in Fig. 6 lies probably in the softness of the central cavity obstacle and also the softness of the effective particle interactions due to the widening of the slit. Apart from that the overall picture of Ref. [26] is confirmed well. In particular, this implies an important observable consequence: The cavity circumference (of approximate length 2πr 0 ) induces a line contribution 2πr 0 γ 2d (ρ 2d,s ) to the insertion free energy of the solute near the wall. Since the insertion free energy is equivalent to the depletion potential up to an additive constant, the latter acquires a term ∝ R(σ − h), the corresponding term in the depletion force is ∝ R/(σ − h). (Note that for h → σ the cavity radius should stay finite as reflected in Fig. 6 (right panel), thus the divergence of the depletion force there according to the geometrical argument of equating the cavity radius with r 0 is unphysical.) The interpretation of this line energy term within the more general framework of morphometric thermodynamics will be given below.
We have seen that in radial direction, the packing of the solvent spheres close to the apex of the annular wedge is well understood by the quasi-2d picture developed above. radial direction. According to Eq. (2), perfect scaling would imply that the depletion force at h = σ is proportional to R. The DFT data violate scaling forr < 0.2 (i.e. when the slit width is smaller than 0.02 σ), as can be seen from Fig. 7 (left panel, inset). This leads to a slow increase of the scaled depletion force f w /(2R), see Fig. 7 (right panel). In contrast, the simulation results of Ref. [33] indicate that f w /(2R) stays constant for large R, quite in agreement with an upper bound derived within the annular slit approximation of Ref. [26] which reads f w /(2R) = −2πγ(ρ s )(1−η 2d ) where η 2d = (π/4) ρ 2d,s σ 2 is the 2d "bulk" packing fraction. Note that the value of the Derjaguin approximation at h = σ, f w /(2R) = −2πγ(ρ s ) is completely off the simulation as well as the DFT values (see next subsection).

Geometric interpretation of the depletion force
Morphological (morphometric) thermodynamics [30] provides a powerful interpretation of the depletion interaction, as has been shown very recently in the special case of the solutesolute interaction [29]. Up to a constant, the depletion potential is nothing but the solvation free energy of the wall and the solute. In morphological thermodynamics the solvation free energy F sol of a body is separated into geometric measures defined by its surface. As the surface, we take the solvent-accessible surface, i.e. the body surface of the combined wallsolute object is defined by the dashed lines in Fig. 2 (a). These geometric measures are the enclosed volume V , the surface area A, the integrated mean and Gaussian curvatures C and X, respectively. To each measure there is an associated thermodynamic coefficient: the pressure p, the planar wall surface tension γ and two bending rigidities κ andκ such that Due to the separation of the solvation free energy into geometrical measures and geometry independent thermodynamic coefficients, it is possible to obtain the coefficients p, γ, κ and κ in simple geometries. The pressure p is a bulk quantity of the fluid. The surface tension γ accounts for the free energy cost of forming an inhomogeneous density distribution close to a planar wall. If the wall is curved the additional free energy cost is measured by κ andκ.
It is possible to obtain all four coefficients from a set of solvation free energies of a spherical particle with varying radius. For a hard-sphere solute in a hard-sphere solvent very accurate analytic expressions for the thermodynamic coefficients are known from FMT-WBII [6], see also App. A.
For the particular case considered here, the depletion potential W (h) follows by subtracting the sum of the solvation free energies of a single wall and a single solute sphere. Thus we find where ∆V and ∆A are the volume and surface area of the overlap of the exclusion (depletion) zones around wall and solute, respectively, and ∆C is the integrated mean curvature of that overlap volume. Note that the fourth characteristic, the integrated Gaussian curvature is the Euler characteristic which is 4π for wall and solute at h ≤ σ (one connected body) and 8π for wall and solute at h > σ (two disconnected bodies). This explains the last term in Eq. (16). The geometric measures ∆V , ∆A and ∆C are given by The depletion force follows as f w (h) = −∂W/∂h. In the limit R → ∞ it is given by: The contribution from the integrated mean curvature is subleading in 1/R but it is quantitatively important since it induces a singularity as h → σ which is proportional to 1/(σ −h) 1/2 . This contribution precisely corresponds to the line tension term associated with the circular apex (of length 2πr 0 ) of the annular wedge identified in the quasi-2d analysis of Ref. [26] and introduced in the last subsection. The line tension in the quasi-2d picture was shown to be γ 2d (ρ 2d,s ) and in the geometric analysis it is −πκ/2. There is good agreement between both expressions [29]. Strictly, the singularity associated with this line term is unphysical of course. Indeed, it can be shown that the morphometric analysis becomes ambiguous very close to h = σ. One would require that the morphometric result for the force is unchanged if another, slightly displaced surface around the solute sphere is chosen for the geometric analysis. In Ref. [29], such an analysis was carried out for the sphere-sphere geometry and the independence of the force on the chosen surface was demonstrated. However, if one chooses e.g. the molecular surface (the surface of the set of points which is never covered by a small solvent sphere), then this surface becomes self-overlapping close to h = σ and one could not expect to ascribe physical significance to surface tensions and mean curvature coefficients of such overlapping surfaces. For the sphere-wall geometry a similar analysis holds. Note that the singularity at h = σ corresponds to a vanishing circumference of the apex (r 0 = 0), however, the 2d analysis of the density distributions at this point yielded a small but finite value for the "effective" apex radius r 0 (see Fig. 6 (right panel)).
In Fig. 8 (left panel) we show FMT-WBII and IE-MV results for the depletion force f w between a big sphere and a wall (size ratio 20, solvent bulk density ρ * s = 0.764) and compare them to MC results from Refs. [31,32] as well as to the morphometric analysis.
The force according to morphometry is plotted only until self-intersection of the molecular surface starts. The FMT data are described very well by the morphometric results; the sharp drop in the depletion force close to h = σ appears to mimick the behavior of the singular term ∝ ∂∆C/∂h. The MC data suffer from relatively large error bars (except for the point h = σ but are overall consistent with both the FMT data and the morphometric analysis. The scatter in the IE-MV data is due to the comparatively low resolution in z-direction as compared to the DFT data (∆z = 0.05 σ vs. ∆z = 0.002 σ). As before, all considered IE closures give similar results, but these correspond to an overall more attractive force compared to DFT and MC. Two other approximations for the depletion force are shown in Fig. 8 (left panel): the DFT insertion route and the Derjaguin approximation. Both approximations fail to describe the depletion force near h = σ. In the insertion route to DFT (see Ref. [20] for its derivation), one circumvents the explicit calculation of the density profile around the fixed wall and solute by exploiting the relation: Here, µ ex b (ρ s ) is the excess chemical potential for inserting one solute into the solvent with density ρ s . The functional derivative of the mixture functional has to be evaluated in the dilute limit for the solutes, ρ b → 0, and with the solvent density profile given by the profile at a hard wall (ρ w (z)). This method relies on an accurate representation of mixture effects for large size ratios in the free energy functional. Therefore, the observed deviations in the depletion force presumably follow from the insufficient representation of higher-order direct correlation functions of strongly asymmetric mixtures in the present forms of FMT [57]. The Derjaguin approximation, on the other hand, corresponds to the leading term for R → ∞ in Eq. (18) inside the depletion region (h < σ), but is easily extended to h ≥ σ [25,26]: with γ slit (h) denoting the excess surface energy of two parallel hard walls at distance h forming a slit. According to Eq. (2), the Derjaguin force scales with R ("colloidal limit") which is an excellent approximation outside the depletion region. (Incidentally, outside the depletion region all methods and the Derjaguin approximation agree with each other.) Inside the depletion region, the neglected mean curvature (or line tension) term is very significant.
Explicit FMT-WBII results for the depletion force for size ratios up to 100 are shown in Fig. 8   expect that for large size ratios, the morphometric picture reliably describes the depletion force. Since the geometry of the overlap volume between two spheres is different from that between a sphere and a wall, we find also a slightly different form for the morphometric force in the limit R → ∞: Thereby, one sees that the scaling relation f w = 2f , valid within the Derjaguin approximation, is violated through the appearance of the mean curvature (line tension) term. For moderate solvent densities up to ρ * s = 0.7 the morphometric form (21) gives a good description of the depletion force for size ratios α 10. This has been reported in Ref. [29]. An analysis of the full solvent density profiles between the two solutes reveals the same features as described in Sec. III A.
It is interesting to investigate the regime where the morphometric description is expected to fail. This will happen when there are strong correlations in the whole annular wedge, i.e. where packing in both the r || -and the z-direction becomes important. Typically, intermediate size ratios and higher solvent densities will induce these strong correlations.
Therefore, we have investigated the depletion force for size ratios 5 and 10 and solvent densities ρ * s = 0.8 and ρ * s = 0.9. In Figs. 9 and 10 we show results for the depletion force and potential at ρ * s = 0.8 for size ratios 5 and 10 respectively. For size ratio 5, the explicit DFT calculations yield a depletion force which does not correspond to the morphometric result at all. It is closer to the insertion route (see (Eq. 19)), though systematically lower.
This yields a depletion potential (Fig. 9 (right panel)) which is about 1 k B T or 25% more attractive at contact than that of the insertion route. For size ratio 10, the explicit DFT data for the depletion force reflect the morphometric form as h → σ, i.e. they give again evidence for importance of the mean curvature (line tension) term in Eq. (21). However, away from that regime, the depletion force deviates significantly from both the morphometric form and the insertion route result such that the depletion potential at contact (Fig. 10 (right   panel)) is about 3 k B T or 50% more attractive at contact than that of the insertion route (for FMT-WBII).
Note that at such a high solvent density (ρ * s = 0.8) there are already significant deviations between FMT-RF and FMT-WBII (Fig. 10 (right panel), circles and squares). Since FMT-WBII has been designed to improve thermodynamic and morphometric consistency at higher densities, the corresponding results can be assumed to be more trustworthy. As an additional check of the numerics, we calculated the depletion force in two ways: (a) via the surface integral in Eq. (1) and (b) via the centered difference of the results for the depletion potential (see Figs. 9 and 10 (left panel), squares and circles). The latter is simply obtained as where the grand potential Ω and the equilibrium density ρ eq (r || , z; h) around the two solutes at distance h are determined by the basic DFT equations (3) and (5) The equality of both routes checks a sum rule check similar to the hard wall sum rule for the problem of one wall immersed in solvent. Weighted-density DFT's usually fulfill the latter [58].
The differences between the insertion route and the explicit FMT data becomes more dramatic for higher densities. In Fig. 11  for ρ * s = 0.9, and the potential well at contact is almost twice as deep as compared with the insertion route. For size ratio 10, the barrier remains (albeit at a different location) but the well depth is equally enhanced by a factor 2 compared with the insertion route as for size ratio 5.
The significant deviations of the depletion potential (according to the explicit FMT-WBII results) from the previously known approximations (Derjaguin approximation and insertion route) sheds some doubts on the quantitative accuracy of simulations of binary hard sphere mixtures (with asymmetries between 5 and 10), employing an effective onecomponent approach for the solutes interacting with their depletion potential. In Ref. [18] the phase diagram of the mixture was scanned in that way using a virial expansion to third order in the Derjaguin approximation for the depletion potential (see the dot-dashed curves in Figs. 9 and 10 (right panel)). Especially for size ratio 5, the explicit DFT results predict a depletion potential with much smaller barrier and a deeper well at contact. This might affect the phase diagram in the corner where the packing fraction of the solutes is low and the solvent packing fraction is high. In Ref. [59] gellation in hard sphere-like, asymmetric mixtures is investigated with an integral equation approach (reference functional approximation [11]) to the depletion potential, akin to the insertion route. We expect quantitative deviations at higher solvent densities compared to explicit DFT calculations which again affect the corner of the phase diagram with low solute and high solvent packing fraction.
Finally we want to mention that the explicit DFT methods introduced here, together with the morphometric form of the depletion force in Eq. (21), can be used to improve available analytic forms [60,61] for the contact values of the pair correlation functions in hard-sphere mixtures (note that the solute-solute contact value is given by g bb (2R b ) = exp[−βW (0)]).

IV. SUMMARY AND CONCLUSIONS
We have presented a detailed analysis of the relationship between depletion forces (solutewall and solute-solute) in a hard solvent and the associated density distributions around the fixed wall/solute particles, using mainly density functional theory of fundamental measure type complemented with Monte-Carlo simulations and integral equation techniques.
For large solute-solvent size ratios, the depletion force is strongly linked to the quasi twodimensional confinement between the solutes (or the solute and the wall). We have shown that the properties of this quasi two-dimensional confined system are reflected in the depletion force by a line-tension term which in turn can be obtained quite generally through morphological thermodynamics. Thus we can formulate an appropriate "colloidal" limit for the depletion force in hard-sphere mixtures: Eq. (18) for the solute-wall case and Eq. (21) for the solute-solute case. This new formulation improves significantly over the Derjaguin approximation which is a frequently employed tool to estimate colloidal interactions in many circumstances.
Although our analysis has been restricted to hard spheres only, the formulation of morphological thermodynamics is very general such that it can be expected that the colloidal limit for the depletion force holds also in mixtures with more general interparticle interactions. However, more explicit tests in this direction are necessary.
For intermediate size ratios 5 and 10 and higher solvent densities ρ * s ≥ 0.8 we found strong deviations in the solute-solute depletion force between the explicit density functional results on the one hand and the Derjaguin approximation or morphometry on the other hand. This is a result of the strong solvent correlations in the annular wedge between the solutes. The depletion well at solute contact is significantly more attractive (by almost a factor of two) than previously estimated. This might have consequences for the phase diagram of binary hard spheres at low solute and high solvent packing fractions. With recently developed techniques [62,63], this question appears to be tractable in direct simulations.
Here, Φ, is a free energy density which is a function of a set of weighted densities {n(r)} = {n 0 , n 1 , n 2 , n 3 , n 1 , n 2 } with four scalar and two vector densities. These are related to the density profile ρ(r) by and the hereby introduced weight functions, depend on the hard sphere radii R i = {R s , R b } of the solvent and solute particles as follows: The excess free energy functional in Eq. (A1) is completed upon specification of the functions ϕ 1 (n 3 ) and ϕ 2 (n 3 ). With the choice we obtain the original Rosenfeld functional [3]. Upon setting we obtain the White Bear functional [4,5], consistent with the quasi-exact Carnahan-Starling equation of state. Finally, with the recently introduced White Bear II functional is recovered.
Next we briefly recapitulate the determination of the thermodynamic coefficients p (pressure), γ (hard wall surface tension), κ andκ (bending rigidity of the integrated mean and Gaussian curvature, respectively) [6]. Consider the free energy of insertion F b for a big solute particle of radius R b and associated radius R = R b + R s of the exclusion sphere around it.
According to Eq. (15), it is given by: On the other hand, F b is equivalent to the excess chemical potential of the solutes in a mixture with solvent particles at infinite dilution (ρ b → 0). This excess chemical potential is obtained as the density derivative ∂/∂ρ b of the mixture free energy (Eq. (A1)), therefore we obtain: By equating the two expressions for F b in eqs. (A7) and (A8), one finds the explicit expressions for the thermodynamic coefficients as linear combinations of the partial derivatives ∂Φ/∂n i . For the Rosenfeld functional, these coefficients are given by: For the White Bear II functional, the thermodynamic coefficients have been given already in Ref. [6] and are reproduced here for completeness: The coefficients of the White Bear II functional are remarkably consistent with respect to an explicit minimization of the functional around solutes of varying radius and fitting the corresponding insertion free energy to Eq. (A7) [6].

APPENDIX B: MINIMIZATION OF FUNDAMENTAL MEASURE FUNCTION-ALS IN CYLINDRICAL COORDINATES
In the following we consider a one-component fluid of solvent hard spheres only (ρ s (r) ≡ ρ(r)). The equilibrium density profile ρ eq (r) of the hard sphere fluid with chemical potential µ = β −1 ln(ρ s Λ 3 ) + µ ex (corresponding to the bulk density ρ s ) in the presence of an arbitrary external potential V (r) is found by minimizing the grand potential which leads to The functional µ[ρ(r)] is given by For the physical problems of this work, sphere-sphere and wall-sphere geometry, the external potential V ext possesses rotational symmetry around the z-axis. Therefore we work in cylindrical coordinates r = (r || cos φ, r || sin φ, z), in which the external potential and the density profile depend only on r || and z, V ≡ V (r || , z) and ρ eq ≡ ρ eq (r || , z). Eqs.

Weighted densities
The convolution integral, defining the weighted densities in Eq. (A2), is calculated by: The 3d Fourier transform of the density profile ρ(r || , z) appearing in the above equation is = FT HT 0 ρ(r || , z) .
Here, we have introduced shorthand notations FT for the Fourier transform in the zcoordinate and HT i for the Hankel transform in the r || -coordinate involving the kernel J i (q || r || ). A vector in real space is given by r = r || e || + ze z , whereas a vector in Fourier space is given by q = q || e ′ || + q z e z , with e || · e ′ || = cos φ q . The 3d Fourier transformsw of the set of weight functions w reduce to sine transforms due to radial symmetry and are explicitly given byw (The remaining weighted densities differ only by a multiplicative factor.) It is convenient to introduce the following parallel and perpendicular components of the Fourier transformed vector weights (k = 1, 2):w Using these definitions, the scalar weighted densities are given by The vector weighted densities, on the other hand, are given by two components (k = 1, 2): n k (r) = n k,|| (r || , z)e || + n k,z (r || , z)e z , We demonstrate this result for the example of the weighted density n 2 . Upon choosing r = (r || , 0, z) we find which is equivalent to Eq. (B16) for k = 2. Here we made use of the integrals Next we consider the evaluation of the convolution type integrals appearing in µ[ρ] (see Eq. (B4)): where p α = ∂Φ/∂n α and α runs over scalar and vector indices. For scalar indices, w k (r) = w k (−r): we recover the standard convolution integral and thus: In the case of vector indices, we observe that the free energy density Φ only depends on vector densities through n 1 · n 2 and n 2 · n 2 . Thus we find (k = 1, 2): p k (r) = p k,|| (r || , z)e || + p k,z (r || , z)e z → (B21) withp k,|| (q || , q z ) = i FT HT 1 p k,|| (r || , z) ,p k,z (q || , q z ) = FT HT 0 p k,z (r || , z) .
Therefore the vector part summands of µ[ρ] are evaluated by We define new variables x = ln(r || /r 0 ) and y = ln(q || /q 0 ) with arbitrary constants r 0 and q 0 .
In terms of these variables whereĴ µ (x) = J µ (q 0 r 0 e x ) andf (x) = e 2x f (e x ). Thus, the Hankel transform takes the appearance of a cross-correlation integral and can be solved via Fourier transforms: The Fourier transform ofĴ µ can be done analytically while the remaining ones are computed numerically using Fast Fourier techniques.
In the actual implementation we followed Ref. [64] which details the proof of orthonormality for discrete functions, defined on a finite interval in logarithmic space (and continued periodically). The following remarks apply: • The fast Hankel transform cannot be applied directly to the density profile, f (r || ) ≡ ρ(r || , z), since it does not go to zero for r || → ∞. It is therefore advantageous to split the external potential into a (possibly) z-dependent background part V (z) and the remainder. The corre- The full equilibrium profile can then be written as ρ eq (r || , z) = ρ V (z)(h(r || , z) + 1) with h(r || , z) → 0 for r || → ∞. The determining equation for h reads ln(h(r || , z) Thus one needs to perform Hankel transforms only on functions which properly go to zero for r || → ∞. In the wall-sphere case, V (z) is naturally given by the wall potential, V s becomes the solute-solvent pair potential u bs and h ≡ h bs is the solute-solvent pair correlation function. In that form, Eq. (B30) resembles the general closure for integral equations (Eq. (10)). In the sphere-sphere case, V (z) = 0 with ρ V = ρ s .
• The assumed periodicity off(x) leads to restrictions on the product q 0 r 0 (low-ringing condition in Ref. [64]). This condition cannot be fulfilled on one grid for both HT 0 (·) and HT 1 (·). However, this does not lead to any noticeable instabilities in the repeated application of the fast Hankel transform.
• In order to avoid aliasing and the amplification of numerical "noise" in the low-x and high-x tails off (x) in the repeated application of the fast Hankel transform we worked with cutoffs r min /r 0 = q min /q 0 = 0.01 and r max /r 0 = q max /q 0 = O(100). The fast Hankel transform itself was calculated on an extended grid x ∈ [−N || ∆x/2, N || ∆x /2] with either N || = 2048, ∆x = 0.01 or N || = 4096, ∆x = 0.005. Outside the "physical" domain defined by x phys ∈ [ln(r min /r 0 ), ln(r max /r 0 )] the functionf (x) was put to zero (a similar prescription applies to y phys ).
Spacing in z-direction was equidistant with grid width ∆z = 0.002 . . . 0.005σ with up to N z = 15, 000 points. Spacing in r || -direction was logarithmic (see above) with N || ≈ In cylindrical coordinates the inhomogeneous OZ equations (7)-(8) read h ss (z, z 0 , s) − c ss (z, z 0 , s) = ds ′ dz ′ ρ V (z ′ )h ss (z, z ′ , |s − s ′ |)c ss (z ′ , z 0 , s ′ ) , h bs (z, z 0 , s) − c bs (z, z 0 , s) = ds ′ dz ′ ρ V (z ′ )h bs (z, z ′ , |s − s ′ |)c ss (z ′ , z 0 , s ′ ) , where the total and direct correlation functions h ij (z, z 0 , s) and c ij (z, z 0 , s) of two particles depend on the distances z and z 0 from the wall and on the projection of their position vectors on the direction along the wall, s = r || − r ||,0 . The integral over s ′ is a 2D convolution, which is most efficiently done by means of the Fast Hankel transform of the zeroth-order (see Appendix B 2). In Fourier space, the OZ equations becomẽ whereh[c](z, z 0 , q || ) = HT 0 h[c](z, z 0 , r || ). The remaining z-integral is evaluated with a simple trapezoidal rule on a uniformly discretized grid of N z points, yielding the following matrix equation: where H and C are N z × N z matrices generated by the corresponding correlation functions at each q || point and R is a diagonal matrix corresponding to ρ V (z)δ(z − z ′ ).
The Lovett-Mou-Buff-Wertheim (LMBW) equation (9) for the background density profile in the presence of a hard wall takes the following form: Here the 2D integration over r || can be considered as the DC component q || = 0 of the Hankel transform of the direct correlation functionc ss (z, z ′ , q || ) and the remaining integral in the second term is taken again with the trapezoidal rule. Note, that the contact density at the plannar wall is known exactly from the wall theorem ρ V (0) = βp, where p is the bulk pressure [66]. Thus, a substitution of the contact density by the expression for the bulk pressure provides the modified LMBW equation for the density profile [67]. Following the ansatz of Plischke and Henderson [68], we used the Carnahan-Starling fit for pressure, which is the first equation in (A10). This in general should improve the accuracy of the density profiles near a plannar wall, if an approximation to the solvent direct correlation function c ss (z, z ′ , r || ) is applied.
The closure equation (10) is the central one in the elaboration of a successful theory, since the OZ and LMBW equations described above are formally exact and need to be complemented with a third relation between between correlation total and direct functions.
Unfortunately, in its exact form this third relation, expressed via additional bridge function b ij (r, r 0 ) is defined only as an infinite sum of highly connected bridge diagrams and so cannot be fully utilized. In practice, one needs to resort to different approximations, which are suited for particular applications depending on the system state. We considered the Percus-Yevick The contact values from the MV closure proved to be in good agreement with the values calculated from Carnahan-Starling equation of state, which constitutes a considerable improvement over the PY and HNC approximations.
In the application of these closures to the inhomogeneous system the following precautions must be observed to obtain robust results.
• Technically we use the same iteration scheme as in Sec. B 3 for the numerical solution of Eqs. (C4)-(C5) and (10) (though due to the memory restrictions the z-grid spacing was only ∆z = 0.05σ with up to N z = 400 points). The convergence drastically depends on the initial conditions, i.e. the zeroth iteration. While solving the PY closure it was enough to initialize both solvent-solvent and solute-solvent correlation functions with zeros, the RY and MV solute-solvent correlations have to be initialized with a good guess to the final solution (which was actually the corresponding PY result). The LMBW equation imposes even stronger requirement on the initial profile, which needs to be either close to the true profile or evolve very slowly from the flat profile together with the corresponding correlation functions.
• It turned out that the order of iterations is crucial for the overall convergency: for a given approximation to the density profile one needs to solve the closure and OZ equations prior to the next iteration on the LMBW equation. It is therefore very costly to get a fully self-consistent solution of all three equations and normally the convergence goal for the background density ρ V (z) is much lower than for the correlation functions.
• To avoid the necessity of requiring the bulk equilibrium behavior in correlation functions for distances far enough from the hard wall, we introduced the second wall at L ∼ 20σ apart from the first one. The wide slit geometry allowed us to keep the zresolution at the maximal possible, yet feasible level. Finer resolution in z in general improves numerical stability of the iteration procedure. Alternatively, one can use a very recent method of the expansion into the orthogonal set of functions proposed by Lado in Ref. [69].