An excess chemical potential for binary hard-sphere mixtures from integral equation theory

In statistical thermodynamics a direct calculation of chemical potential is a difficult task than computing other quantities, such as an internal energy, the pressure, and the correlation functions determining the structures of the physical systems. Various theoretical methods have been extensively employed to calculate the chemical potential, such as, the Widom’s test particle insertion method Widom63 or its modifications Adams74 ; Attard93 and the scaled-particle Monte-Carlo (SP-MC) method Barosova96 ; Borowko00 ; Martin18 .

In theoretical studies for liquid an integral equation (IE) theory has been widely used due to its low computational cost and a fact that the theory is a thermodynamically complete, when combined with an appropriate closure equation Ornstein14 ; Percus58 ; Chandler72 ; Hirata81 ; Perkyns92 ; Kovalenko99 ; Luchko16 ; Lebowitz58 ; Tsogbayar19 . From the IE approach, calculating the chemical potential is almost straightforward employing an analytical expression based on the correlation functions. Therefore, in this work our purpose is to calculate the chemical potential for binary hard-sphere mixtures using the IE approach with the Percus-Yevick closure [8]. We will also compare our numerical findings with the SP-MC results obtained by Barošová et. al., Barosova96 .

The paper consists of four sections. The first section presents the Introduction, while the second and the third sections discuss a brief description of theoretical formulation and numerical results, respectively. Then a conclusion follows.

For the atomic mixture system, the site-site Ornstein-Zernike equation (SSOZ) Chandler72 ; Hirata81 ; Perkyns92 ; Kovalenko99 which establishes a relation between the total correlation function $h_{ij}$ and the direct correlation function $c_{ij}$ in the $\mathbf{k}$-space can be written in the form

where $\rho_{i}=x_{i}\rho$ is a number density, and $x_{i}$ is the mole fraction of component $i$, and $\rho=\sum_{i}\rho_{i}$ is the total number density. The tildes in Eq. (1) denote a Fourier transformation. Since we do not know both correlation functions in Eq. (1), we cannot solve this equation directly. Therefore, to solve this equation, we need an auxiliary equation which we call a closure relation, and then they can be solved self-consistently. A general closure for the mixture may be given in the form

Here $u_{ij}$ is an interaction potential, which in this work is given by

where $\sigma_{ij}=\frac{1}{2}(\sigma_{i}+\sigma_{j})$, and $\sigma_{i}$ is a diameter of component $i$, and $\beta=1/k_{\mathrm{B}}T$, $k_{\mathrm{B}}$ is a Bolztmann’s constant and $T$ is the temperature for the system. .

The Percus-Yevick (PY) bridge function Percus58 which we will use in our calculation has a form

where $\gamma_{ij}=h_{ij}-c_{ij}$ is the indirect correlation function. We note that in this PY approximation a solution for an IE for the hard-sphere systems can be found in an analytical form Lebowitz58 . Once we have the correlation functions describing the structure of the mixture, we can compute an excess chemical potential for component $i$ using a following expression

Note that a derivation of this expression can be performed similarly as it was obtained for the pure Lennard-Jones fluid by Tsogbayar and Luchko Tsogbayar19 . Moreover, the expression can be used for any bridge functions since it does not require the explicit forms of them.

In our calculation we chose a component 1 as a reference particle which is a larger component of the mixture, that is $(\sigma_{1}>\sigma_{2})$, and a packing fraction $\eta$ is given as

where $x_{2}=1-x_{1}$ is the mole fraction for component 2. We solve an Eq. (1) numerically using the Picard iteration method. We performed numerical calculations for four different values of diameter ratio $\sigma_{2}/\sigma_{1}:0.3,0.5,0.6,0.9$ and at $\eta=0.49$. This value of the packing fraction indicates that our calculation can be done in the stable liquid region since the value is slightly lower that the freezing packing fraction for pure hard spheres Barosova96 . We do our calculation as function of the mole fraction $x_{1}$ of the larger component. In all calculations a number of grid points is $N=2^{15}$, and a length interval is $L=32\sigma_{1}$. The convergence criterion for a successive two iterations of the indirect correlation function is set to $10^{-6}\sigma_{1}$.

The columns 2 and 3 in Table I show our numerical results of an excess chemical potential $\beta\mu^{ex}_{1}$, for various values of the mole fraction $x_{1}$ and at $\eta=0.49$. We placed results of the SP-MC method Barosova96 in the columns 3 and 5. When $\sigma_{2}/\sigma_{1}=0.3$, and for the smaller values of $x_{1}$ our calculated values show small discrepancies, however they reduce with an increase of the mole fraction of the larger component. For $\sigma_{2}/\sigma_{1}=0.9$, our values placed in column 4 are close to the SP-MC data shown in the last column. We note that the number of particles used in the conventional canonical ensemble Monte-Carlo method of Ref. Barosova96 is $1728$, however, since our IE method is an implicit method, it is independent on number of particles involved in the system, that is, we do not care about number of particles in the IE formulation.

In Table II we present a comparison of our calculated values for $\beta\mu^{ex}_{2}$ with available values obtained by the SP-MC method Barosova96 for $\sigma_{2}/\sigma_{1}=0.3$ and $0.9$, and at $\eta=0.49$. All our findings in this table are close to an accurate SP-MC values Barosova96 .

Once we compute the excess chemical potential, we may obtain directly the logarithm of the activity coefficient using a following equation Barosova96 :

In Fig. 1 we show a logarithm of the activity coefficient calculated with an Eq. 6 for the four different diameter ratios as function of $x_{1}$ at constant $\eta=0.49$ for the smaller component. Blue, red, green and black curves correspond to the diameter ratios of 0.3, 0.5 0.6 and 0.9, respectively. Looking at these curves, the logarithm of the activity coefficient decreases very quickly with $x_{1}$, especially for the smaller diameter ratios in which the spheres’ diameters are more asymmetric, and for $\sigma_{2}/\sigma_{1}=0.9$, its feature looks almost a linear (black). The crosses are the SP-MC values taken from Ref. Barosova96 . The curves and crosses follow each other.

Fig. 2 presents same plots as shown in Fig. 1, but for the larger component. For this case, the logarithm of the activity coefficient is negative. It decreases relatively quickly for smaller diameter ratios as the mole fraction of the larger component increases, which is indeed a similar pattern appeared in Fig. 1. Crosses showing the SP-MC values Barosova96 follow our obtained curves.

Red curves $(\sigma_{2}/\sigma_{1}=0.5)$ in Fig. 1 and Fig. 2 are not compared with any data. However, looking at other obtained curves, we can note that they can be as accurate as others.

In this work we have solved the site-site Ornstein-Zernike equation for binary hard-sphere mixture for four different diameter ratios and at $\eta=0.49$ in the PY approximation. The excess chemical potential for two components as a function of the mole fraction $x_{1}$ of the larger component has been computed with an analytical expression containing the correlation and bridge functions, which in general can be used in an evaluation of the excess chemical potential from the IE approach with any forms of bridge function. Moreover, the logarithms of the activity coefficients for the components as a function of the mole fraction $x_{1}$ have been obtained and compared with accurate SP-MC data. It has been shown that our obtained results agree well with those data.