Recent developments in classical density functional theory: Internal energy functional and diagrammatic structure of fundamental measure theory

An overview of several recent developments in density functional theory for classical inhomogeneous liquids is given. We show how Levy's constrained search method can be used to derive the variational principle that underlies density functional theory. An advantage of the method is that the Helmholtz free energy as a functional of a trial one-body density is given as an explicit expression, without reference to an external potential as is the case in the standard Mermin-Evans proof by reductio ad absurdum. We show how to generalize the approach in order to express the internal energy as a functional of the one-body density distribution and of the local entropy distribution. Here the local chemical potential and the bulk temperature play the role of Lagrange multipliers in the Euler-Lagrange equations for minimiziation of the functional. As an explicit approximation for the free-energy functional for hard sphere mixtures, the diagrammatic structure of Rosenfeld's fundamental measure density unctional is laid out. Recent extensions, based on the Kierlik-Rosinberg scalar weight functions, to binary and ternary non-additive hard sphere mixtures are described.


Introduction
The theoretical study of inhomogeneous classical liquids received a boost through the development of classical density functional theory (DFT). Evans' 1979 article [1] constitutes a central reference to DFT. The proof of the variational principle that underlies the theory, including the existence and the uniqueness of the free energy functional, can be viewed as the classical analogue of Mermin's earlier (1965) work on quantum systems at finite temperatures [2]. This forms a generalization of the Hohenberg-Kohn theorem for ground state properties of quantum systems. Historic milestones that predate these developments, and that can be re-formulated in classical DFT language, are Onsager's 1949 treatment of the isotropicnematic liquid crystal phase transition in systems of long and thin hard rods [3], and van der Waals' 1893 theory of the microscopic structure of the liquid-gas interface [4].
Both the Hohenberg-Kohn proof and the Mermin-Evans proof start from a variational principle for the respective many-body function. In the quantum case, this is the Rayleigh-Ritz inequality for the manybody groundstate wave function. In the classical case, the theory rests on the Gibbs inequality for the equilibrium many-body phase space distribution. The corresponding functionals are the ground state energy in the quantum case and the thermodynamic grand potential in the classical case. Both functionals depend (trivially) on the position-dependent external one-body potential. Via an intricate sequence of arguments [1,2], which has become textbook knowledge [5], the dependence on the external potential is played back to the more useful dependence on the (in general) position-dependent one-body density distribution.
In 1979 Levy showed that quantum DFT can be obtained in a more compact and straightforward way via a method that he called the constrained search [6] (see also [7]). Here, the search for the minimum is performed in the space of all trial many-body wave functions. The constraint that is fixed during this search is that all trial wave functions considered need to generate the given one-body density distribution. Levy's derivation is both rigorous and elegant and constitutes a standard reference for electronic structure DFT. Among the impressive number of citations of his 1979 article, there are only very few papers that draw connections to classical DFT, [8] being an example, although the method permits rather straightforward application to the classical case [9]. See also Percus' general concept of "overcomplete" density functionals [10].
The constrained search method offers two significant benefits over the Mermin-Evans proof. One is simplicity, avoiding the reductio ad absurdum chain of arguments. The other is that it yields an explicit definition of the Helmholtz free energy density as a functional of the (trial) one-body density distribution. This formula, as reproduced in (2.2) below, is explicitly independent of the external potential. The generalization that facilitates this development is the concept of minimization in the space of all phase space distributions under the constraint of a given one-body density (2.4). The issue of representability of the one-body density has been addressed in [11]. It turns out that the simplicity of Levy's method allows one to construct more general DFTs. An example is the variational framework developed in [12], which rests on the internal energy functional (rather than the Helmholtz free-energy functional), which depends on the one-body density and on a local (position-dependent) entropy distribution. A dynamical version of this theory, based on linear irreversible thermodynamics and phenomenological reasoning, is proposed in [12]. Several common approximations for free energy functionals were transformed to internal energy functionals. Having reliable approximations for the functional is a prerequisite for applying DFT to realistic (three-dimensional) model systems. The task of constructing usable functionals is different from the conceptual work outlined so far. In particular, a calculation of the constrained search expressions for the (free-energy or internal-energy) functional would amount to an exact solution of the many-body problem. Hence, the importance of these expressions is rather of conceptual nature.
In 1989 Rosenfeld wrote a remarkable letter in which he proposed an approximate free energy functional for additive hard sphere mixtures [13]. His theory unified several earlier liquid state theories, such as the Percus-Yevick integral equation theory for the bulk structure [5], the scaled-particle theory for thermodynamics, and Rosenfeld's own concepts, such as the scaled-field particle theory of [13], and encapsulated these into what he called fundamental measure theory (FMT). Kierlik and Rosinberg in 1990 [14] re-wrote the same functional [15] in an alternative way, using only four (not six, as Rosenfeld) weight functions to build weighted densities via convolution with the bare one-body density of each hard sphere species. The two strands of FMT were pursued both with significant rigour and effort, see the recent reviews [16][17][18]. Rosenfeld's more geometric approach was extended to further weight functions by Tarazona [19] in his treatment of freezing. A critical discussion of the properties of the relevant convolution kernels can be found in [20]. Very recently, Korden [21,22] demonstrated the relationship of the FMT with the exact virial expansion.
It is vexing that Rosenfeld himself, who was certainly very versed in the diagrammatic techniques of liquid state theory (see e.g. [23]), apparently neither analysed nor formulated his very own FMT within such a framework. The non-local structure of FMT, its coupling of space integrals via convolution, and the plentiful appearances of the one-body density distribution(s) seem to constitute an ideal playground for formulation in diagrammatic language. A comprehensive understanding of the diagrammatic nature of FMT could not help just to ascertain and clarify the nature of the approximations that are involved [21], but also, and more importantly from a pragmatic point of view, could enable one to construct new functionals for further model systems. Taking the FMT weight functions as bonds in a diagrammatic formulation, one immediately faces the combinatorial problems associated with their number, which in the Rosenfeld-Tarazona formulation are at least seven per hard sphere species (four scalar, two vector, one tensor), which makes the book-keeping task a seemingly daunting one.
It was shown [24,25] that one can formulate the Kierlik-Rosinberg form of FMT in a diagrammatic way. The concept was applied to one-dimensional hard rods, where it gives Percus' exact result [26,27], as well as to five-dimensional hard hypersphere mixtures, where it gives a functional that outperforms Percus-Yevick theory, and improves on previous FMT attempts [28]. Two crucial properties of the diagrammatic formulation can be identified. One is the relative simplicity of the diagrams that describe the coupling of the various space integrals in the density functional. The topology of the diagrams is of star-like or tree-like shape, hence providing a significant reduction as compared to the complexity of the exact virial series. The number of arms is equal to the power in density and the bonds are weight functions rather than Mayer functions, like they are in the exact virial expansion. The book-keeping problem of having to deal with a large number of different weight functions is addressed and rendered almost trivial by exploiting the tensorial structure that underlies the Kierlik-Rosinberg form of FMT [29], where the geometric index of the weight functions is a proper tensorial index, and corresponding (isometric and metamorphic) transformations can be formulated [29]. Hence, the fully scalar Kierlik-Rosinberg formulation turned out to be indeed simpler to handle, and easier to generalize for full control of the degree of non-locality in the functional. These developments facilitated the generalization of FMT for binary non-additive hard spheres [30] to ternary mixtures [31].
In the present contribution, we describe the basic ideas underlying the above developments, without the full detail that is given in the respective original papers, but with further illustrative examples in order to provide an introduction to the subject. The paper is organized as follows. In section 2, Levy's method is sketched, and both the free-energy and the internal-energy functionals are defined. A basic introduction to the diagrammatic formulation of FMT is given in section 3, including a brief overview of applications to non-additive hard sphere fluids in bulk and at interfaces. Section 4 gives concluding remarks.

Variational principle for the grand potential functional
Consider a classical many-body system with position coordinates r 1 , r 2 , . . . , r N and momenta p 1 , p 2 , . . . , p N , where N is the number of particles. Denote the classical trace in the grand ensemble by . . dr N dp 1 . . . dp N , where h is Planck's constant. For a given total interatomic potential U (r 1 , . . . , r N ) between the particles, one defines the intrinsic Helmholtz free energy functional by the explicit expression [9] F where T is temperature, k B is the Boltzmann constant, m is the particle mass, p i = |p i |, and the minimization searches all many-body (phase space) probability distributions f (r 1 , . . . , r N , p 1 , . . . , p N ; N ) that are normalized according to Tr f = 1, (2.3) and that yield the fixed trial one-body density ρ(r) via The relationship (2.4) is indicated as f → ρ in the notation of (2.2). It can be shown [9] that the intrinsic free energy functional (2.2) satisfies the Euler-Lagrange equation (2.5) where µ is the chemical potential, v(r) is an external one-body potential acting on the system, so that the total external potential for a given microstate is N i=1 v(r i ), and ρ 0 (r) is the equilibrium one-body density,

43603-3
where the equilibrium many-body distribution f 0 is given by the Boltzmann distribution, Here, the normalization constant is the grand partition sum and the Hamiltonian is A concise version of the proof of (2.5) is laid out in the following.

Sketch of the constrained search proof a la Levy
It can easily be shown [1] that the equilibrium grand potential Ω 0 for Hamiltonian H N is obtained from minimizing Mermin's form of the grand potential functional, where f is a trial many-body distribution. Decompose the right hand side into a double minimization where the inner minimization is a search under the constraint that the f generates ρ via (2.4). For suitable Hamiltonians as (2.9) with external potential v(r), this can be written as In the expression above where F [ρ] is given by (2.2) and hence in equilibrium the Euler-Lagrange equation (2.5) follows. Finally, note that from (2.11) the grand potential density functional can be defined as

DFT based on the internal-energy functional
Using two constraints, rather than one, one can define the internal-energy functional as +U (r 1 , . . . , r N ) , (2.16) where the constraint for the density is (2.4) and that for the local entropy distribution s(r) is The equilibrium value for the local entropy is with the equilibrium many-body distribution f 0 given in (2.7), and the total entropy obtained as the spatial integral S 0 = drs 0 (r). The Euler-Lagrange equations have the form Explicit forms of various internal-energy functionals can be found in [12], as can be a dynamic prescription, similar in spirit to dynamical DFT [32], for the joint time evolution of s(r, t ), ρ(r, t ), where t is time.
The proof of (2.16)-(2.20) is a straightforward application of Levy's method, as laid out in section 2.2; for the explicit form see [12].

Diagrammatic formulation for additive and non-additive mixtures
We turn to hard sphere mixtures and use the common splitting of the intrinsic Helmholtz free energy into ideal gas and excess (over ideal gas) contributions where ρ i (r) is the one-body density distribution, Λ i is the thermal de Broglie wavelength of species i , and F exc [{ρ i }] is the Helmholtz excess free energy functional that is due to interparticle interactions. Its lowest order (in density) virial expansion is where the Mayer function is defined as being the pair interaction potential between species i and j . The sum in (3.2) is over all species i , j . For hard spheres f i j (r < σ i j ) = −1 and zero otherwise; here σ i j is the hard core interaction distance between species i and j . Kierlik and Rosinberg introduced a set of four "weight functions" w 0 (R i , r ), w 1 (R i , r ), w 2 (R i , r ), and w 3 (R i , r ), where R i = σ ii /2 is the particle radius, and r is radial distance [14]. Using the weight functions, the hard sphere Mayer function can be written as a sum of convolution integrals,

43603-5
where the asterisk denotes the convolution of two functions h 1 (r ) and h 2 (r ), which is defined as (h 1 * h 2 )(|r − r ′ |) = dxh 1 (r − x)h 2 (r ′ − x). The spatial arguments of the weight functions have been omitted in the notation of the right hand side of (3.3). The direct space expressions of the weight functions are is the Dirac distribution, and Θ(·) is the Heaviside step function.
The appearance of particular combinations of products in (3.3) can be based on dimensional analysis. In order to see this, note that the weight functions w ν are objects with dimensions of (length) ν−3 , hence w 0 , w 1 , w 2 , and w 3 carry dimensions of (length) −3 , (length) −2 , (length) −1 , and (length) 0 , respectively. Each of the products in (3.3) has the dimension (length) −3 , which cancels the (length) 3 which is due to the convolution integral and hence yields a dimensionless Mayer function.
This property can be formalized [33] in order to establish more mathematical structure. One can re-write (3.3) as (3.4) where the coefficients M ντ are chosen in such a way that only the "allowed" combinations of weight functions contribute. Hence, for most of the index combinations ντ, the coefficients vanish, M ντ = 0. The crucial step is to view the coefficients M ντ as the elements of a metric M which operates in the space of weight functions w(R,r). Here, a vector of weight functions w is defined by its components w ν (R, r), hence w = (w 0 , w 1 , w 2 , w 3 ). The length scale R acts as a parameter. As an aside, one should not confuse the vector w with Rosenfeld's vectorial weight functions w v1 (r) and w v2 (r) [13], which have a very different origin, rooted in the geometry of the sphere in three-dimensional space. The w(R,r), on the other hand, are elements of a four-dimensional vector space; their index runs from 0 to 3.
The metric M can be represented as a matrix, which is akin to a mirrored 4 × 4 unit matrix in that its   (3.6) where the superscript t indicates matrix transposition, and the spatial arguments have been omitted in the notation; these are the same as in (3.4).
In order to illustrate the framework, we display corresponding diagrams for binary hard sphere mixtures in figure 1. The deconvolution of the Mayer bond f 11 (r ) between particles of species 1 is displayed in figure 1 (a). The length of the w bonds indicates the magnitude of the argument R, in this case R 1 . The kink is located at position x, which is the integration variable in the convolution integral (3.4). Hence, in this and in the following diagrams, the position of a kink (or junction to be introduced below) is integrated over. The open circles indicate fixed positions r and r ′ , i.e., the arguments in f 11 (|r − r ′ |). One commonly refers to these as root points. Recall that multiplying each root point by the one-body density distribution and integrating over its position yields (up to a factor of k B T /2) the i j = 11 contribution to the exact low density limit of the excess free energy functional (3.2). Figure 1 (b) gives the corresponding diagram for species 2, here taken to be of a larger size (hence R 2 > R 1 ). The cross species diagram, i j = 12, is shown in figure 1 (c). From the diagrammatic representation it is clear that the total length of the diagram, i.e., the range over which f 12 (r ) is non-zero, is fully determined by the accumulated length of the two arms, R 1 and R 2 . Hence, the cross species interaction distance is σ 12 = R 1 + R 2 , a case to which one refers to as an additive hard sphere mixture. Additivity is a fundamental property of the space of weight functions, and is indeed more general than outlined so far. Consider a given vector of weight functions w(R,r), where the length scale R is fixed. We aim at finding a linear operation K(d), where d is a length scale, that yields

43603-6
where * represents the combined operation of spatial convolution ( * ) and dot product (·) in the fourvector space. Hence, in a less concise but more explicit way (3.7) can be written componentwise as (3.8) where the components of K are denoted by K τ ν and the two indices ν and τ run from 0 to 3. That such a K exists is a non-trivial matter, and forms the heart of the binary non-additive hard sphere functional [30], see the detailed investigation of the properties of K in [33]. The existence and properties are most easily demonstrated in Fourier space representation, where the convolutions become mere products. The Fourier transform of the "kernel matrix" K possesses a representation as the matrix exponential K(R, q) = exp (RG) , (3.9) where the "generator" is (3.10) As an almost trivial consequence of (3.9), the relationship

43603-7
holds, simply due to the exponential satisfying exp((R + R ′ )G) = exp(RG) exp(R ′ G). Correspondingly, in real space Two features of this expression are welcome. One is that the the Kierlik-Rosinberg weight functions w ν (R, r ) constitute four of the sixteen (4 × 4) components of K(r,R). Hence, the Kierlik-Rosinberg weight functions are "automatically" generated when starting with (3.9) and (3.10). The second feature is that all further components K τ ν (R, r ) share their degree of non-locality with the weight functions, i.e., all components of K vanish for distances r > R. This becomes a crucial property when one uses these objects as more general convolution kernels in the construction of third and higher orders (in density) in the excess free energy functional. Before we describe such developments we return to the problem of a general binary hard sphere mixtures, in which the cross interaction is non-additive, i.e., σ 12 R 1 + R 2 = (σ 11 + σ 22 )/2. Equation (3.12) allows us to use two convolutions rather than only one in order to represent the cross Mayer function for a non-additive binary mixture as − f 12 (r ) = w(R 1 ) * K(d) * w t (R 2 ). (3.13) Here, the length scale d is due to the non-additivity, and satisfies σ 12 = R 1 + d + R 2 . Figure 1 (d) shows the corresponding diagrams, where K(d) constitutes an additional type of bond. The length of the bond can be adjusted, via changing d, in order to generate the (given) range of f 12 (r ). In the diagram, the position of each kink is integrated over. These integrals correspond to the convolution integrals in (3.13). The low density expansion of the binary non-additive hard sphere functional [30] contains the same interspecies diagrams as the additive Kierlik-Rosinberg functional, as shown in figure 1 (a) and (b), and contains the diagram shown in figure 1 (d) for the cross species interaction.
Using the property (3.12), one can go beyond (3.13) and represent the single K matrix via further deconvolution as a (convolution) product of two K matrices, with appropriately chosen length scales d 1 and d 2 , so that σ 12 = R 1 + d 1 + d 2 + R 2 , and hence − f 12 (r ) = w(R 1 ) * K(d 1 ) * K(d 2 ) * w t (R 2 ). (3.14) The corresponding diagram is shown in figure 1 (e). This constitutes the cross species low-density limit of the ternary non-additive hard sphere functional [31]. For ternary mixtures, the three cross species lengthscales σ 12 , σ 13 , σ 23 are decomposed as (3.15) with appropriate values of d 1 , d 2 and d 3 , which can be uniquely determined for a mixture with three components [31].
We turn to the third-virial level. The exact contribution to the excess free energy functional is Due to the pairwise coupling of the three integration variables, one commonly refers to (3.16) as (the sum of) triangle diagrams. FMT fails to generate this exact expression [20], but yields a very reasonable approximation to it. Apart from going from Mayer bonds as convolution kernels to weight function bonds, the crucial step is a "re-wiring" or topological change [21] in the structure of the diagrams. Avoiding the loop in the triangle diagram (a problem that gets more severe with an increasing order in density in the virial expansion, see below), the FMT diagrams connect the (three) weight function bonds to a common central space integral. A tree-like (or star-like) topology results, see figure 2. Writing out the space integrals explicitly yields . (3.17) From an algebraic point of view, the operation on the right hand side of (3.17) is similar to a triple scalar product of w(R i ), w(R j ), and w(R k ). The third-rank tensor J that generates this operation can be derived [24] from the sole requirement that the result, after the x-integral in (3.17), does not contain unphysical divergences (more precisely that the corresponding order in the partial bulk direct correlation function is finite everywhere). The tensor J is symmetric in all its indices; its non-vanishing components are J 033 = J 123 = 1, and J 222 = 1/(4π).  K is a second rank tensor, which is contracted via an M metric to a w bond and, on its other side to the third-rank junction tensor J. This allows one to control the range of the cross Mayer bond f 12 (r ) freely. The w bonds are specific for each species (indicated by the species index 1,2 at each root point).

43603-9
Figure 2 illustrates these diagrams for the case of an additive binary mixture. Note that, as on the second virial level, all bond lengths are pre-determined, once the particle radii R 1 and R 2 have been chosen. These are uniquely determined from R i = σ ii /2, leaving no room to go beyond the additive case. See [20] for an in-depth discussion of the so-called "lost cases" of FMT that come from the intrinsic differences between the triangle and the three-arm star diagrams. In short, when stretching all bonds maximally, the exact diagram is larger than the FMT approximation; see [20]. Nevertheless, the different lengths of the FMT arms approximate the equilateral cases i j k = 111, figure 2 (a) and i j k = 222, figure 2 (d), as well as simultaneously the isosceles cases i j k = 112, figure 2 (b) and i j k = 122, figure 2 (c). For non-additive mixtures, however, this is is insufficient for approximating all triangle diagrams.
For non-additive mixtures, K can be used in order to replace one of the weight functions w(R) in (3.17) by K(d) * w(R). The resulting diagrams are shown in figure 3. The intra-species contributions, i j k = 111 in figure 3 (a) and i j k = 222 in figure 3 (d), are unchanged as compared to the additive case shown in figure 2. The inter-species diagrams differ by an additional K bond that acts as a spacer between the arm of the "minority"-species (i.e., the one that appears only once, not twice) and the central space integral.
Note that both "majority" arms directly connect to the center. This applies to i j k = 112, as shown in figure 3 (b), and to i j k = 122 as shown in figure 3 (c).
Although the FMT expressions for third-order diagrams are approximations, they possess one important property which is most clearly analysed when considering the low-density limit of the (partial) bulk fluid two-body direct correlation functions, c i j (r ), obtained as second functional derivatives,  (3.19) where the contribution linear in densities is obtained from the triangle diagram by multiplying one of the root points with the bulk density and integrating its position over space, i.e., turning the root point of species k into a density field point. The field is constant in the homogeneous bulk. A selection of the corresponding integrals is displayed in figure 4. The FMT results give the exact result for c * i j k (r ) for all combinations of species. This applies to all versions of the theory, whether additive [13,14], binary non-additive [30], or ternary non-additive [31], and is a first indication of successful description of bulk structure. The construction of the additive hard sphere FMT free energy functional can be based on a diagrammatic series using the assumption that all diagrams of higher than third order in density are also of star-like shape [24]. The weight functions w that constitute the arms are connected to a single central space integral. In order to connect p arms of the p-th order in density, a tensor of rank p is required.
One can show [24] that a recursive relation holds and that these object can be obtained by tensor contraction of an appropriate number (p − 2) of third-rank tensors J. As an illustration, at fourth order where the first tensor on the right hand side has one index lowered via contraction with the metric, Together with the coefficients 1/[p(p − 1)], which are taken from the zero-dimensional properties of the system [24], the resulting series is displayed in figure 5. The series can be explicitly summed, and yields, in three spatial dimensions, the Kierlik-Rosinberg form of FMT. Each arm together with its filled field point constitutes a "weighted density" [13,14], . (3.20) The diagrams that constitute the FMTs for binary and ternary non-additive hard sphere mixture, possess different topology; a summary is displayed in figure 6. The one-component hard sphere FMT functional, shown in figure 6 (a) consists of stars with arms of the same length R. In one dimension (hard rods on a line), this is equivalent to Percus' exact functional [26]. The binary version possesses two different types of arms corresponding to the small (here species 1) and large (species 2) component, see figure 6 (b). The one-dimensional version of the functional for non-additive hard core mixtures [34] is 43603-10  The bonds are fundamental measure weight functions that are joined by an M metric (second order in density), or by junction tensor J of k-th rank (k-th order in density). The position of the inner junction (center of the star) is integrated over. This corresponds to the "outer" integral dx over the free energy density Φ(n τ (x)) in the standard, weighted density formulation with weighted densities n τ (x). The scalar coefficient of the k-th order diagram is 1/[k(k − 1)].
an approximation that differs from the exact solution [35]. The binary non-additive FMT has two central junctions, where the arms of each of them belong to one species only (the small species 1 belongs to the left center, and the large species 2 to the right center). The two centers are connected via a K bond of appropriate length d, so that σ 12 = R 1 + d + R 2 , see figure 6 (c). Finally, figure 6 (d) displays a typical diagram of the ternary non-additive FMT. Here, there are three subcenters, one for each species, with corresponding arms. The subcenters are connected via K bonds to one global center. The positions of all junctions are integrated over. A full account of the ternary FMT functional can be found in [31].

Non-additive hard sphere fluids in bulk and at interfaces
We give a brief summary of applications of the FMT for binary non-additive hard sphere mixtures. Based on the FMT for this system [30,31] but also on Monte Carlo computer simulations, a range of physical phenomena were investigated. The bulk fluid structure on the two-body level is described very satisfactorily by the theory, as compared to simulation data. This applies to the results for the partial pair correlation functions, g i j (r ), when obtained through the Ornstein-Zernike route using the partial direct correlations functions, c i j (r ), obtained as the second functional derivative (3.18) of the free energy functional. This is a quite severe test of the theory, because both the derivative and the Ornstein-Zernike equations constitute involved operations. The general performance of the theory is very satisfactory, although the core condition g i j (r < σ i j ) = 0 is only approximately satisfied. This can be rectified with Percus' test-particle method [36], where the grand potential is minimized in the presence of an external potential that is equal to the inter-particle interaction. Test-particle results for the g i j (r ) reproduce the simulation data very well [37], as long as the system is away from the "decoupling case" of negative nonadditivity, so that σ 12 ≪ (σ 11 + σ 22 )/2. It is important to note that the theory does not yield unphysical artifacts in test-particle results for partial pair correlation functions [37,38] The theory yields analytic expressions for the partial direct correlation functions (both in real space and in Fourier space), and hence, via the Ornstein-Zernike relation also analytic expressions for the partial structure factors S i j (q). This makes it very convenient to carry out an analysis of the asymptotic, large distance decay of bulk pair correlation functions from pole analysis of the S i j (q) in the plane of complex wavevectors q [39,40], and to relate these results to the decay of one-body density profiles in inhomogeneous situations, such as at interfaces [41].

43603-11
For positive non-additivity, σ 12 > (σ 11 + σ 22 )/2, which is sufficiently large, the system displays fluidfluid phase separation into two fluid phases with different chemical compositions. The theory gives good account of the location of the fluid-fluid demixing binodal in the plane of partial packing fractions of the two species [30,41]. A wealth of interesting interfacial phenomena results as a consequence of the bulk fluid demixing, and we refer the reader to the original papers on fluid demixing, asymptotic decay of correlations and free fluid interfaces [41], first-order layering and critical wetting transitions in nonadditive hard sphere mixtures [42], and capillary condensation of non additive hard sphere mixtures in planar confinement [43].

Conclusions
We have described a range of recent developments and applications of classical density functional theory in the study of bulk and interfacial properties of liquids. Starting with a reassessment of the underlying variational principle, we have laid out how to use Levy's constrained search method in order to define the free-energy functional. Levy's method provides us with an explicit expression for the freeenergy functional (2.2). We showed that the concept can be generalized in order to define an internalenergy functional, which possesses the one-body density distribution and a local entropy distribution as trial fields. A dynamical theory built on the internal-energy functional can be found in [12]. Using Levy's method, the definition of the functional is explicitly independent of the external potential, which consti-tutes a conceptual advantage. However, the minimization in the function space of many-body probability distributions cannot be in practice carried out for a realistic model Hamiltonian. Hence, one has to rely on approximations for the functional, as is common in DFT. We refer the reader to [12] for a description of various approximate internal-energy functionals, which have been obtained from Legendre transforming the corresponding approximation for the Helmholtz free-energy functional.
We have laid out the basic ideas underlying the recent progress in formalizing the mathematical structure of FMT. This includes the tensorial nature, in that the "geometric" indices of the Kierlik-Rosinberg weight functions play the role of tensorial indices. A diagrammatic notation helps to clarify the non-local nature of the excess free energy functional. We have given an overview of applications of the binary non-additive hard sphere functional to a range of phenomena.
The above work extends the previous efforts that were primarily based on the intimate connection of the properties of the free energy functional under dimensional crossover, i.e., the result of applying the functional to density distributions that correspond to an extreme confinement in one or more spatial directions, made it possible to generalize FMT to the Asakura-Oosawa-Vrij model [44,45] of colloid-polymer mixtures [46], the Widom-Rowlinson model [47], and penetrable spheres that interact with a constant repulsive plateau [48]. These models have in common that their zero-dimensional properties, i.e., the statistical mechanics of a cavity of the size of a single particle, carries still some of the essentials of the true three-dimensional problem. Hence, the zero-dimensional problem, which can be solved exactly (or with good approximations [48]) in the above cases, is sufficient as a central modification over the FMT for hard spheres, in order to obtain a reliable free energy functional. See e.g., [49][50][51][52][53] for applications to confined model colloid-polymer mixtures.
In future work it would be interesting to apply the dynamical test particle limit [54,55] to non-additive hard sphere mixtures in order to gain a better understanding of the dynamical behaviour of such mixtures. A further important problem is to re-consider formulating an FMT for soft sphere models [56][57][58][59] in the light of the diagrammatic and tensorial structure. Work along these lines is in progress [60]. It would also be interesting to investigate the implications of Levy's method for the statistical mechanics of quenched-annealed fluid mixtures, where DFT was obtained both via the replica trick [61][62][63][64], and via a first-principles derivation following the Mermin-Evans arguments [65].