Chemical-potential route for multicomponent fluids

It is well known that in a fluid made of particles interacting via a pair-wise potential, the most relevant statistical-mechanical quantity in equilibrium is the pair correlation function, or radial distribution function (RDF), $g(r)$ Barker and Henderson (1976); Hansen and McDonald (2006). It is defined as the average number density at a distance $r$ from a certain reference particle, relative to the global density. Apart from accounting for the structural correlation properties of the fluid, the RDF $g(r)$ allows one to obtain a number of thermodynamic quantities $\{\psi_{i}\}$ by means of elegant formulas. Since those thermodynamic quantities are connected by differential relations, one would expect to get the same macroscopic description, i.e., the same free energy $A$, regardless of the specific route $g(r)\to\psi_{i}\to A$ followed not . However, this is not necessarily the case when an approximate RDF is used as a starting point, what results in the well-known thermodynamic inconsistency problem Barker and Henderson (1976).

To be more specific, let us assume an $s$-component system made of $N=\sum_{\alpha=1}^{s}N_{\alpha}$ particles ($N_{\alpha}$ being the number of particles of species $\alpha$) enclosed in a $d$-dimensional volume $V$. The total number density is $\rho=N/V$ and the partial number densities are $\rho_{\alpha}=N_{\alpha}/V=x_{\alpha}\rho$, $x_{\alpha}=N_{\alpha}/N$ being the mole fractions. The pair interaction potential, not necessarily isotropic, between two particles of species $\alpha$ and $\gamma$ is $\phi_{\alpha\gamma}(\mathbf{r})=\phi_{\gamma\alpha}(-\mathbf{r})$. In terms of the pair correlation function of species $\alpha$ and $\gamma$, $g_{\alpha\gamma}(\mathbf{r})$, it is possible to express the pressure $p$ as Hansen and McDonald (2006)

where $Z$ is the compressibility factor and $\beta\equiv 1/k_{B}T$ ($k_{B}$ and $T$ being the Boltzmann constant and the temperature, respectively). Analogously, the internal energy per particle $u=U/N$ can be expressed as

Equations (1) and (2) are the virial (or pressure) and energy routes to thermodynamics, respectively. The free energy $A=Na$ (where $a$ is the free energy per particle) can be (partially) obtained not either from Eq. (1) or from Eq. (2) by inversion of the thermodynamic relations

Apart from the virial and energy routes, a popular route is the compressibility one. It reads

where the element $\widehat{h}_{\alpha\gamma}$ of the matrix $\widehat{\mathsf{h}}$ is proportional to the zero wavenumber limit of the Fourier transform of the total correlation function $h_{\alpha\gamma}(\mathbf{r})=g_{\alpha\gamma}(\mathbf{r})-1$, namely

Combination of Eqs. (3) and (5) allows one to get the free energy from the RDF via the compressibility route.

Equations (3) and (4) are related to the derivatives of the free energy per particle $a$ with respect to the total number density and to the temperature, respectively, but they ignore the composition-dependence of $a$. In fact, the partial derivatives of $a$ with respect to the partial number densities are related to the chemical potentials:

While in Eq. (3) the free energy per particle $a$ is seen as a function of the total number density $\rho$ and the $s-1$ independent mole fractions $\{x_{1},x_{2},\ldots,x_{s-1}\}$, in Eq. (7) $a$ is seen as a function of the $s$ partial densities $\{\rho_{1},\rho_{2},\ldots,\rho_{s}\}$. Adopting the former point of view, Eq. (7) can be rewritten as

This allows us to get the useful identity

Recognizing that $\left[\partial(\rho a)/{\partial\rho}\right]_{\beta,\{x_{\alpha}\}}=a+p/\rho$ is the Gibbs free energy per particle, Eq. (9) is not but the fundamental equation of thermodynamics Callen (1960).

Cross partial derivatives of the free energy should be independent of the order of derivation, which yields the following Maxwell relations involving the chemical potential:

In Eqs. (11) and (12) the replacement $\rho_{\alpha}\to N_{\alpha}/V$ must be done in order to express the chemical potential as a function of $\{N_{\alpha}\}$ and $V$.

The interesting question is, how can one obtain the chemical potential $\mu_{\nu}$ (and hence the free energy) from the knowledge of the RDF? This question has been addressed in several textbooks Hill (1956); Reichl (1980) and classical papers Mandell and Reiss (1975); Reiss et al. (1959), but the most general formula (valid for any dimensionality, interaction potential, and coupling protocol) seems not to have been derived yet. The first aim of this paper is the derivation of that chemical-potential route to thermodynamics [see Eq. (36) below]. The second goal is the application of the route to a mixture of (additive) hard spheres (HS) in the context of the Percus–Yevick (PY) approximation [see Eq. (64) below]. This extends to multicomponent systems, the work recently reported in Ref. Santos (2012).

This paper is organized as follows. The formal derivation of the chemical-potential route is carried out in Sec. II. As a simple test, it is checked in Sec. III that the use of the exact pair correlation function to first order in density provides the exact second and third virial coefficients. Section IV is devoted to the specialization of the chemical-potential route to $d$-dimensional HS mixtures with zero or positive nonadditivity. In the three-dimensional case, the knowledge of the contact values of the radial distribution function for additive HS mixtures in the scaled-particle theory (SPT) and PY approximations is exploited to obtain the chemical potential. It is observed that, while the SPT yields a result thermodynamically consistent with the virial route, this is not the case of the PY result. Finally, the paper ends with some concluding remarks in Sec. VI.

As stated before, we consider an $s$-component mixture with $N_{\alpha}$ ($\alpha=1,\ldots,s$) particles of species $\alpha$ and $N=\sum_{\alpha=1}^{s}N_{\alpha}$ total number of particles in a volume $V$. We will employ the short-hand notations $\mathbf{N}\equiv\{N_{1},\ldots,N_{s}\}$, $\mathbf{r}^{N}\equiv\{\mathbf{r}_{1},\ldots,\mathbf{r}_{N}\}$, and $\text{d}\mathbf{r}^{N}\equiv\text{d}\mathbf{r}_{1}\cdots\text{d}\mathbf{r}_{N}$, $\mathbf{r}_{i}$ being the spatial coordinates of particle $i$. If we denote by $\Phi_{\mathbf{N}}(\mathbf{r}^{N})$ the total potential energy, the (canonical-ensemble) configurational probability density is Hill (1956); Hansen and McDonald (2006)

where

is the configurational integral. The average number of pairs of particles (per unit volume) of species $\alpha$ and $\gamma$ located at $\mathbf{r}_{\alpha}$ and $\mathbf{r}_{\gamma}$, respectively, is

where $\epsilon_{i}$ denotes the species particle $i$ belongs to. We define the pair correlation function as $g_{\alpha\gamma}(\mathbf{r}_{\alpha},\mathbf{r}_{\gamma})=n_{\alpha\gamma}(\mathbf{r}_{\alpha},\mathbf{r}_{\gamma})/\rho_{\alpha}\rho_{\gamma}$. Inserting Eq. (14) into Eq. (16), one obtains

where, without loss of generality, particles $i=1$ and $j=2$ are assumed to belong to species $\alpha$ and $\gamma$, respectively.

Contact with thermodynamics is made through the free energy $A_{\mathbf{N}}$ of the system:

where the partition function $\mathcal{Z}_{\mathbf{N}}$ factorizes as

being the ideal-gas partition function. In Eq. (20), $\Lambda_{\alpha}=h/\sqrt{2\pi m_{\alpha}k_{B}T}$ (where $h$ is the Planck constant and $m_{\alpha}$ is the mass of a particle of species $\alpha$) is the thermal de Broglie wavelength.

Let us now focus on a given species $\nu$. The associated chemical potential is

where

In the second line of Eq. (23) we have taken into account that $N_{\nu}\gg 1$ and it is understood that the subscript $\mathbf{N}+1$ means $\{N_{1},\ldots,N_{\nu}+1,\ldots,N_{s}\}$. Obviously, $Q_{\mathbf{N}+1}(\beta,V)$ is given by Eq. (15) with the replacements $N\to N+1$ and $\mathbf{N}\to\mathbf{N}+1$. Moreover, without loss of generality, we will assign the label $i=0$ to the extra particle of species $\nu$, so that $\text{d}\mathbf{r}^{N+1}=\text{d}\mathbf{r}_{0}d\mathbf{r}^{N}$.

Now we assume that the potential energy is pair-wise additive, namely

where $\phi_{\alpha\gamma}(\mathbf{r}_{\alpha},\mathbf{r}_{\gamma})=\phi_{\alpha\gamma}(\mathbf{r}_{\gamma}-\mathbf{r}_{\alpha})$ is the interaction potential of two particles of species $\alpha$ and $\gamma$.

In order to establish a relationship between the chemical potential and the pair correlation functions, it is convenient to introduce a coupling parameter $\xi$ such that its value $0\leq\xi\leq 1$ controls the strength of the interaction of particle $i=0$ to the rest of particles. A similar charging process was employed by Onsager in a different context Onsager (1933). It is also closely related to the so-called Widom insertion method Widom (1963).

In our specific problem, the interaction potential between particles $i=0$ and $j\geq 1$ is $\phi_{\nu\epsilon_{j}}^{(\xi)}(\mathbf{r}_{0},\mathbf{r}_{j})$, with the boundary conditions

The associated total potential energy and configuration integral are

Note that Eqs. (27) and (28) (with $0<\xi<1$) define a system of $N+1$ particles and $s+1$ species, one of the species being made by the single particle $i=0$. Only in the limits $\xi=0$ and $\xi=1$ does one recover a number of species $s$. Analogously to Eq. (17), the pair correlation function of particle $i=0$ and any particle of species $\alpha$ is

Obviously,

Assuming that the potentials $\phi_{\nu\alpha}^{(\xi)}$ are differentiable with respect to $\xi$, the second line of Eq. (23) can be rewritten as

Now, from Eqs. (27) and (28) we obtain

where again particle $j=1$ is assumed (without loss of generality) to belong to species $\alpha$. Making use of Eq. (29),

Combination of Eqs. (33) and (35) yields the desired result

where we have taken into account the translational invariance property $\phi_{\nu\alpha}^{(\xi)}(\mathbf{r}_{\nu},\mathbf{r}_{\alpha})=\phi_{\nu\alpha}^{(\xi)}(\mathbf{r}_{\alpha}-\mathbf{r}_{\nu})$. In terms of the cavity function $y_{\nu\alpha}^{(\xi)}(\mathbf{r})=g_{\nu\alpha}^{(\xi)}(\mathbf{r})e^{\beta\phi_{\nu\alpha}^{(\xi)}(\mathbf{r})}$, Eq. (36) becomes

Equation (36) or, equivalently, Eq. (37) constitutes the chemical-potential route to thermodynamics. It expresses the chemical potential of species $\nu$ in an $s$-component mixture (the “solvent”) in terms of the pair correlations of an $(s+1)$-component mixture where the extra component consists of a single particle (the “solute”) whose interaction with the rest of the particles is controlled by a coupling parameter $\xi$, which switches from $\xi=0$ (the solute does not feel the presence of the solvent particles) to $\xi=1$ (the solute is indistinguishable from a solvent particle of species $\nu$). The result is independent of the number of species $s$, the dimensionality of the system $d$, the pair interaction functions $\phi_{\alpha\gamma}(\mathbf{r})$ (which can be isotropic or not), and the protocol $\phi_{\nu\alpha}^{(\xi)}(\mathbf{r})$ followed to go from $\phi_{\nu\alpha}^{(0)}(\mathbf{r})=0$ to $\phi_{\nu\alpha}^{(1)}(\mathbf{r})=\phi_{\nu\alpha}(\mathbf{r})$. If the exact $g_{\nu\alpha}^{(\xi)}(\mathbf{r})$ is employed, one gets the exact chemical potential $\mu_{\nu}$ (and, hence, the exact free energy $A$) regardless of the protocol. This is illustrated in Sec. III at the level of the third virial coefficient. On the other hand, if an approximate function $g_{\nu\alpha}^{(\xi)}(\mathbf{r})$ is used instead, not only the resulting approximate free energy will be different from the one derived from any of the other routes (e.g., virial, energy, and compressibility), but it will, in general, depend on the choice of the protocol. Thus, the chemical-potential route adds an extra source of thermodynamic inconsistency.

The aim of this section is to exploit Eq. (37) when the pair correlation function is truncated to first order in density. First, let us rewrite Eq. (37) as

where $f_{\nu\alpha}^{(\xi)}(\mathbf{r}_{\nu},\mathbf{r}_{\alpha})\equiv e^{-\beta\phi_{\nu\alpha}^{(\xi)}(\mathbf{r}_{\nu},\mathbf{r}_{\alpha})}-1$ is the Mayer function Hansen and McDonald (2006); Hill (1956).

To first order in density,

where $f_{\gamma\alpha}(\mathbf{r}_{\gamma},\mathbf{r}_{\alpha})\equiv e^{-\beta\phi_{\gamma\alpha}(\mathbf{r}_{\gamma},\mathbf{r}_{\alpha})}-1$. Insertion into Eq. (38) yields

where we have called

Exchanging $\alpha\leftrightarrow\gamma$ and $\mathbf{r}_{\alpha}\leftrightarrow\mathbf{r}_{\gamma}$ in the summation and the integral, $I_{\nu}^{(\xi)}$ can be rewritten as

Therefore, Eq. (40) becomes

where

and the translational invariance property has been applied again. Note that both $\mathcal{F}_{\nu\alpha}$ and $\mathcal{G}_{\nu\alpha\gamma}$ are invariant under any permutation of indices. This guarantees that Eq. (13) is verified.

Equations (43)–(45) show that, as expected, the specific protocol does not play any role in the final result. From Eq. (7) we finally get

Making use of Eq. (3) it is straightforward to get

where the second and third virial coefficients are

As expected, these are the exact expressions Hansen and McDonald (2006).

Let us now particularize the chemical-potential route (37) to HS mixtures. In that case,

where $r=|\mathbf{r}|$, $\Theta(x)$ is Heaviside’s step function, and $\sigma_{\alpha\gamma}$ is the range of the infinitely repulsive interaction between particles of species $\alpha$ and $\gamma$.

We now need to introduce the extra particle $i=0$ (the solute) coupled to the remaining $N$ particles through a coupling parameter $0\leq\xi\leq 1$ via the set of interaction potentials $\phi_{\nu\alpha}^{(\xi)}(\mathbf{r})$. According to Eq. (26), $\exp[-\beta\phi_{\nu\alpha}^{(0)}(\mathbf{r})]=1$ and $\exp[-\beta\phi_{\nu\alpha}^{(1)}(\mathbf{r})]=\Theta(r-\sigma_{\nu\alpha})$, but otherwise a certain freedom to fix the protocol $\phi_{\nu\alpha}^{(\xi)}(\mathbf{r})$ exists. The most natural choice is a HS form, i.e.,

with $\sigma_{\nu\alpha}^{(0)}=0$ and $\sigma_{\nu\alpha}^{(1)}=\sigma_{\nu\alpha}$. Therefore,

so that Eq. (37) becomes

where we have taken into account that the integral of $\text{d}\mathbf{r}$ over all orientations is $d2^{d}v_{d}r^{d-1}\text{d}r$, $v_{d}=(\pi/4)^{d/2}/\Gamma(1+d/2)$ being the volume of a $d$-dimensional sphere of unit diameter. Note that in Eq. (53) the integration variable has changed from $\xi$ to $\sigma_{\nu\alpha}^{(\xi)}$ and we have simplified the notation as $\sigma_{\nu\alpha}^{(\xi)}\to\sigma_{0\alpha}$ and $y_{\nu\alpha}^{(\xi)}\to y_{0\alpha}$. This change avoids the need to specify the $\xi$-dependence of $\sigma_{\nu\alpha}^{(\xi)}$.

Equation (53) applies to any choice of the set $\{\sigma_{\alpha\gamma}\}$. Henceforth we will assume the condition $\sigma_{\alpha\gamma}\geq\frac{1}{2}(\sigma_{\alpha}+\sigma_{\gamma})$, where $\sigma_{\alpha}$ and $\sigma_{\gamma}$ are the diameters of particles of species $\alpha$ and $\gamma$, respectively. This implies that the HS mixture is either additive or has a positive nonadditivity. In such a case, Eq. (53) can be further simplified. To that end, it is convenient to decompose $\beta\mu_{\nu}^{\text{ex}}$ into two pieces,

where

In the first contribution, the condition $\sigma_{0\alpha}\leq\frac{1}{2}\sigma_{\alpha}$ implies that the solute can lie “inside” a particle of species $\alpha$. This makes the contribution $\beta\mu_{\nu}^{\text{ex},\text{I}}$ very easy to evaluate. First, by reversing the steps leading from Eq. (23) to Eq. (53), one can write

where $\xi_{0}$ denotes the common value of the coupling parameter $\xi$ at which $\sigma_{\nu\alpha}^{(\xi)}=\sigma_{0\alpha}$ takes the value $\frac{1}{2}\sigma_{\alpha}$ for all $\alpha$. According to Eqs. (27), (28), and (51),

Next, thanks to the condition $\sigma_{\alpha\gamma}\geq\frac{1}{2}(\sigma_{\alpha}+\sigma_{\gamma})$, we note that in Eq. (58) the spatial regions excluded to particle $i=0$ by the remaining particles $j=1,\ldots,N$ do not overlap, so that

where $\eta=v_{d}\sum_{\alpha}\rho_{\alpha}\sigma_{\alpha}^{d}$ is the total packing fraction of the solvent system. Thus, Eq. (58) reduces to $Q_{\mathbf{N}+1}^{(\xi_{0})}(\beta,V)=Q_{\mathbf{N}}(\beta,V)(1-\eta)$ and therefore $\beta\mu_{\nu}^{\text{ex},\text{I}}=-\ln(1-\eta)$.

Taking all of this into account, we finally get

In this section, we apply Eq. (60) to the derivation of the chemical potential of a three-dimensional additive HS fluid mixture, i.e., $\sigma_{\alpha\gamma}=\frac{1}{2}(\sigma_{\alpha}+\sigma_{\gamma})$, as resulting from the PY and SPT approximations. Given an $s$-component mixture, the corresponding contact values are

where

and the parameter $q$ takes the values $q=0$ and $q=\frac{3}{4}$ for the PY Lebowitz (1964) and SPT Reiss et al. (1959); Helfand et al. (1961); Lebowitz et al. (1965); Mandell and Reiss (1975); Rosenfeld (1988); Heying and Corti (2004) theories, respectively. The more accurate Boublík–Grundke–Henderson–Lee–Levesque (BGHLL) Boublík (1970); Grundke and Henderson (1972); Lee and Levesque (1973) expression corresponds to the intermediate value $q=\frac{1}{2}$.

Now, in order to apply Eq. (60), we assume that a solute particle $i=0$ is introduced in such a way that it interacts with a particle of species $\alpha$ through a HS potential of range $\sigma_{0\alpha}$, the contact value of the corresponding cavity function being $y_{0\alpha}(\sigma_{0\alpha})$. In principle, we are free to choose $\sigma_{0\alpha}$ within the interval $\frac{1}{2}\sigma_{\alpha}\leq\sigma_{0\alpha}\leq\sigma_{\nu\alpha}$. On the other hand, if we want to make use of the approximation (61), we need to make the solute particle interact additively with the rest of the particles in the system. This is achieved by assuming that the solute particle is a sphere of diameter $\sigma_{0}$ within the range $0\leq\sigma_{0}\leq\sigma_{\nu}$ and $\sigma_{0\alpha}=\frac{1}{2}(\sigma_{0}+\sigma_{\alpha})$. In that case, Eq. (60) becomes

where $y_{0\alpha}(\sigma_{0\alpha})$ is simply given by Eq. (61) with $\sigma_{\gamma}\to\sigma_{0}$. After simple algebra, Eq. (63) yields