Quantum Perturbation Theory in Fluid Mixtures

1. Introduction

Experimental assessment of macroscopic thermo-dynamical parameters under extreme conditions is almost impossible and very expensive. Therefore, theoretical EOS for further experiments or evaluation is inevitable. In spite of other efficient methods of calculation such as integral equations and computer simulations, we have used perturbation theory because of its extensive qualities. Moreover, other methods are more time consuming than perturbation theories. When one wants to deal with realistic intermolecular interactions, the problem of deriving the thermodynamic and structural properties of the system becomes rather formidable. Thus, perturbation theories of liquid have been devised since the mid-20th century. Thermodynamic perturbation theory offers a molecular, as opposed to continuum approach to the prediction of fluid thermodynamic properties. Although, perturbation predictions are not expected to rival those of advanced integral-equations or large scale computer simulations methods, they are far more numerically efficient than the latter approaches and often produced comparably accurate results.

Dealing with light species such as $$He$$ and $$H_2$$ at low temperature and high densities makes it necessary taking into account quantum mechanical effects. Quantum rules and shapes related with the electronic orbital change completely the macroscopic properties.

Furthermore, for this fluid mixture, the quantum effect has been exerted in terms of first order quantum mechanical correction term in the Wigner-Kirkwood expansion. This term by generalizing the Wigner-Kirkwood correction for one component fluid to binary mixture produce acceptable results in comparison with simulation and other experimental data. Since utilizing Wigner-Kirkwood expansion in temperatures below 50 K bears diverges, we preferred to restrict our investigations in ranges above those temperatures from 50 to 4000 degrees. In these regions our calculations provide more acceptable results in comparison with other studies.

This term make a negligible contribution under high temperatures conditions. Taking into account various contributions, we have utilized an improved version of the equation of state to study the Helmholtz free energy F, to investigate the effects of P and T on thermodynamic properties of helium and hydrogen isotopes mixtures over a wide range of densities. We also have studied effects of concentrations of each component on macroscopic parameters. In addition, comparisons among various perturbation and ideal parts have been presented in logarithmic diagrams for different densities and concentrations for evaluation of perturbation terms validity in respect to variables ranges.

The first section is dedicated to a brief description of Wigner expansion which leads to derivation of first quantum correction term in free energy. With the intention of describing effects of quantum correction term we have explained theoretical method of our calculations in the frame work of statistical perturbation theory of free energy in section two. In section three we have depicted diagrams resulted from our theoretical evaluations and gave a brief explanation for them. In section four we have focused on the description of our calculations and its usages in different areas. Finally, some applications of this study have been introduced in the last section.

2. Quantum correction term

Considering quantum system of $$N$$ identical particles of mass $$m$$ confined to the region of $$\Lambda$$ with the interacting potential of$$U$$. This structure is considered in $$\upsilon$$-dimension space ($$R^{\upsilon }$$). In the absence of external fields the Hamiltonian of particles is given as

$$H=\frac{1}{2m}{\left(-i\hslash \overrightarrow{\nabla }\right)}^2+U\left(\overrightarrow{r}\right)$$E1

Where, $$\hslash$$ is the Plank constant. The equilibrium statistical mechanics of the particle system is studied in the canonical ensemble at the temperature $$T$$ (or, alternatively, the inverse temperature $$\beta =1/k_{B}T$$ with $$k_B$$ being Boltzmann’s constant). Quantum effects will be considered via de Broglie wavelength $$\lambda =\hslash \sqrt{\beta /m}$$. For a typical microscopic length of particles $$l$$, for sufficiently small dimensionless parameter $$\lambda /l$$ semi-classical regime is dominant. In such system Boltzmann density in configuration space $$\overrightarrow{r}$$ can be expanded in powers of $${\lambda }^2$$ within the well-known Wigner-Kirkwood expansion [1, 2]. In the case of an inverse-power-law repulsive potential $$V\left(r\right)=V_0{\left(a/r\right)}^n$$ from the range $$1<n<\infty$$, the Wigner-Kirkwood expansion turns out to be analytic in $${\lambda }^2$$ [3]. In the hard-core limit $$n\to \infty$$, this expansion is not further correct and one has the non-analyticity of type $${\left({\lambda }^2\right)}^{1/2}$$, as was shown in numerous analytic studies [4-7]. In contrast to the bulk case, the resulting Boltzmann density involves also position dependent terms which are non-analytic in $$\lambda$$. Under some condition about the classical density profile, the analyticity in $$\lambda$$ is restored by integrating the Boltzmann density over configuration space.

2.1. Wigner-Kirkwood expansion

To have an analytical equation for quantum effects in fluid we must derive partition function of it. In approximating partition function we need to evaluate Boltzmann density. Consequently having an expansion of quantum correction terms it is necessary to expand Boltzmann density. Considering system of $$N$$ particles in the infinite space in standard Wigner-Kirkwood expansion [1, 2] fermions or boson exchange effects between quantum particles have been neglected. In the “bulk” regime, equilibrium quantities of this system in the nearly classical regime can be expanded in powers of $$h^2$$. In this section, we review briefly the derivation of this expansion for utilizing it in statistical perturbation framework. The Boltzmann density $$B_{\beta }$$ in configuration space $$\overrightarrow{r}$$ can be formally written in the basis of plane waves as a $$\upsilon N$$-dimensional integral defined in an infinite domain $$R^{\upsilon }$$:

$$B_{\beta }=⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩={\displaystyle \int \frac{d\overrightarrow{P}}{{\left(2\pi \hslash \right)}^{\upsilon N}}{e}^{-\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}{e}^{-\beta H}{e}^{\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}}$$E2

Where $$\overrightarrow{p}=\left(p_1,p_2,p_3,\mathrm{...}\right)$$ is the $$\upsilon N$$-dimensional momentum vector. Instead of considering we take the Laplace transform of this operator with respect to the inverse temperature $$\beta$$,

$${\displaystyle {\int }_0^{\beta }d\beta e^{-\beta H}}{e}^{-\beta z}=\frac{1}{H+Z}$$E3

Via integrating equation 2 in respect to $$\beta$$ we have

$${\displaystyle {\int }_0^{\beta }d\beta ⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩e^{-\beta z}=}{\displaystyle \int \frac{d\overrightarrow{P}}{{\left(2\pi \hslash \right)}^{\upsilon N}}{e}^{-\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}\frac{1}{H+Z}{e}^{\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}}$$E4

Let us introduce following definition

$$H+z=D+Q$$E5

That $$Q$$ and $$D$$ respectively represent

$$\begin{array}{l}Q=\frac{1}{2m}{\left(-i\hslash \overrightarrow{\nabla }\right)}^2-\frac{1}{2m}{\overrightarrow{P}}^2\\ D=\frac{1}{2m}{\overrightarrow{P}}^2+U\left(\overrightarrow{r}\right)+z\end{array}$$E6

One can expand

$$\frac{1}{H+Z}=\frac{1}{D}-\frac{1}{D}Q\frac{1}{D}+\frac{1}{D}Q\frac{1}{D}Q\frac{1}{D}-\mathrm{...}$$E7

$$Q$$, operates in the following manner

$$Q\left[f\left(\overrightarrow{r}\right)e^{\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}\right]=e^{-\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}\left[\frac{i\hslash }{m}\overrightarrow{p}.\overrightarrow{\nabla }+\frac{{\hslash }^2}{2m}{\overrightarrow{\nabla }}^2\right]f\left(\overrightarrow{r}\right)$$E8

And then we can find that

$$\begin{array}{l}{e}^{-\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}\frac{1}{H+z}{e}^{\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}=\frac{1}{D}{\displaystyle \sum _{n=0}^{\infty }{\left\{\left[\frac{i\hslash }{m}\overrightarrow{p}.\overrightarrow{\nabla }+\frac{{\hslash }^2}{2m}{\overrightarrow{\nabla }}^2\right]\frac{1}{D}\right\}}^n}\\ {\displaystyle \int \frac{d\overrightarrow{P}}{{\left(2\pi \hslash \right)}^{\upsilon N}}{e}^{-\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}\frac{1}{H+z}{e}^{\left(i/\hslash \right)\overrightarrow{p}.\overrightarrow{r}}}={\displaystyle \int \frac{d\overrightarrow{P}}{{\left(2\pi \hslash \right)}^{\upsilon N}}\frac{1}{D}{\displaystyle \sum _{n=0}^{\infty }{\left\{\left[\frac{i\hslash }{m}\overrightarrow{p}.\overrightarrow{\nabla }+\frac{{\hslash }^2}{2m}{\overrightarrow{\nabla }}^2\right]\frac{1}{D}\right\}}^n}}=\\ {\displaystyle {\int }_0^{\beta }d\beta ⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩e^{-\beta z}=}{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \int \frac{d\overrightarrow{P}}{{\left(2\pi \hslash \right)}^{\upsilon N}}\frac{1}{D}}{\left\{\left[\frac{i\hslash }{m}\overrightarrow{p}.\overrightarrow{\nabla }+\frac{{\hslash }^2}{2m}{\overrightarrow{\nabla }}^2\right]\frac{1}{D}\right\}}^n}\end{array}$$E9

So we have expanded series in $${\hslash }^{2n}$$ which enable us power series of $${\hslash }^n$$. It remains to define $$1/D^j$$

$$\frac{1}{D^j}={\displaystyle {\int }_0^{\infty }d\beta \frac{1}{\left(j-1\right)!}}{\beta }^{j-1}{e}^{-\beta D}={\displaystyle {\int }_0^{\infty }d\beta e^{-\beta z}\frac{1}{\left(j-1\right)!}}{\beta }^{j-1}{e}^{-\beta \left[{\overrightarrow{p}}^2/2m+U\left(\overrightarrow{r}\right)\right]}$$E10

and finally integrating on the momentum variables $$\overrightarrow{p}$$, the Boltzmann density in configuration space is obtained as the series

$$⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩={\displaystyle \sum _{n=0}^{\infty }{B}_{\beta }^{\left(n\right)}\left(\overrightarrow{r}\right)},$$E11

where

$$\begin{array}{l}{B}_{\beta }^{\left(0\right)}\left(\overrightarrow{r}\right)=\frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U}\\ B_{\beta }^{\left(1\right)}\left(\overrightarrow{r}\right)=\frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U}{\lambda }^2\left[\frac{-\beta }{4}{\overrightarrow{\nabla }}^{2}U+\frac{{\beta }^2}{6}{\left(\overrightarrow{\nabla }U\right)}^2\right],\\ B_{\beta }^{\left(2\right)}\left(\overrightarrow{r}\right)=\frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U}\left\{{\lambda }^2\left[\frac{\beta }{6}{\overrightarrow{\nabla }}^{2}U-\frac{{\beta }^2}{8}{\left(\overrightarrow{\nabla }U\right)}^2\right]+O{\lambda }^4\right\}\\ \text{We conclude that the quantum Boltzmann density in configuration space is given, }\\ \text{to order}{\lambda }^2\text{, by}\\ ⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩=\frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U}\left\{1+{\lambda }^2\left[\frac{-\beta }{12}{\overrightarrow{\nabla }}^{2}U+\frac{{\beta }^2}{24}{\left(\overrightarrow{\nabla }U\right)}^2\right]+O{\lambda }^4\right\}=\\ \frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}\left\{e^{-\beta U}\left[1-\frac{{\lambda }^2\beta }{24}{\overrightarrow{\nabla }}^{2}U\right]+\frac{{\lambda }^2}{24}{\overrightarrow{\nabla }}^2{e}^{-\beta U}+O{\lambda }^4\right\}\end{array}$$E12

Integrating Boltzmann density ignoring exchange effects over configuration space will result in partition function of fluids mixture.

$$Z_{qu}=\frac{1}{N!}{\displaystyle {\int }_{V}d\overrightarrow{r}⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩}$$E13

Substituting the $$\lambda$$-expansion of the Boltzmann density (12A) into formula (13), the quantum partition function takes the expansion form

$$Z_{qu}=\frac{1}{N!}{\displaystyle {\int }_{\Lambda }d\overrightarrow{r}\frac{1}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}\left\{e^{-\beta U}\left[1-\frac{{\lambda }^2\beta }{24}{\overrightarrow{\nabla }}^{2}U\right]1+\frac{{\lambda }^2}{24}{\overrightarrow{\nabla }}^2{e}^{-\beta U}+O{\lambda }^4\right\}}$$E14

For expressing macroscopic physical quantities, one defines the quantum average of a function $$f\left(\overrightarrow{r}\right)$$ as follows

$${⟨f⟩}_{qu}=\frac{1}{Z_{qu}N!}{\displaystyle {\int }_{\Lambda }d\overrightarrow{r}⟨\overrightarrow{r}|e^{-\beta H}|\overrightarrow{r}⟩}f\left(\overrightarrow{r}\right)$$E15

At the one-particle level, one introduces the particle density

$$n_{qu}\left(r\right)={⟨{\displaystyle \sum _{j=1}^N\delta \left(r-r_j\right)}⟩}_{qu}$$E16

At the two-particle level, the two-body density is given by

$$n_{qu}^{\left(2\right)}\left(r,r^{\prime }\right)={⟨{\displaystyle \sum _{\begin{array}{l}j,k=1\\ j\ne k\end{array}}^N\delta \left(r-r_j\right)\delta \left(r^{\prime }-r_j\right)}⟩}_{qu}$$E17

And the pair distribution function

$$g_{qu}\left(r,r^{\prime }\right)=\frac{n_{qu}^{\left(2\right)}\left(r,r^{\prime }\right)}{n_{qu}\left(r\right)n_{qu}\left(r^{\prime }\right)}$$E18

The classical partition function and the classical average of a function $$f\left(\overrightarrow{r}\right)$$ are defined as follows

$$Z=\frac{1}{N!}{\displaystyle {\int }_{\Lambda }\frac{d\overrightarrow{r}}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U\left(\overrightarrow{r}\right)}}$$E19

$$⟨f⟩=\frac{1}{ZN!}{\displaystyle {\int }_{\Lambda }\frac{d\overrightarrow{r}}{{\left(\sqrt{2\pi }\lambda \right)}^{\upsilon N}}{e}^{-\beta U\left(\overrightarrow{r}\right)}f\left(\overrightarrow{r}\right)}$$E20

Consequently with the definition of equation 19 one can derive below equation for $$Z_{qu}$$

$$Z_{qu}=Z\left\{1-{\lambda }^2\frac{\beta }{24}⟨{\overrightarrow{\nabla }}^{2}U⟩+O{\lambda }^4\right\}$$E21

$$\beta F_{qu}=-\ln \left(Z_{qu}\right)$$E22

$$\ln \left(Z_{qu}\right)=\ln (Z)+\ln (1-{\lambda }^2\frac{\beta }{24}⟨{\overrightarrow{\nabla }}^{2}U⟩+O{\lambda }^4+\mathrm{...})$$E23

Since we have $$\ln (1-x)=-x-x^2/2-\mathrm{...}$$ we can expand the second term in the right side. By means of equation 18 in deriving $$⟨{\overrightarrow{\nabla }}^{2}U⟩$$ we can have explicit formula for the second term of which indicates the first term of Wigner-Kirkwood correction part that is consist of the second derivative of potential function that leads to below equation for quantum correction term with the number density of $$n$$ we have

$$F_{qu}^{(1)}=\frac{h^2{N}_{A}n\beta }{96{\pi }^{2}m}{\displaystyle {\int }_{{\sigma }^0}^{\infty }{\nabla }^{2}U(r)g(r)dr}$$E24

$$g(r)$$ represents radial distribution function, which is a measure of the spatial structure of the particles in reference system, is the expected number of particles at a distance $$r$$. $$N_A$$ is Avogadro constant and $${\sigma }^0$$ is the distance in which potential function effectively tend to zero.

2.2. Free energy

Generalizing to multi-component system we have [8]

$$F_{qu}^{(1)}=\frac{h^2{N}_{A}n\beta }{96{\pi }^2}{\displaystyle \sum _{i,j}\frac{c_i{c}_j}{m_{ij}}{\displaystyle {\int }_{{\sigma }_{ij}^0}^{\infty }{\nabla }^2{u}_{ij}(r)g_{ij}(r)4\pi r^2\overline{V_{ij}}dr}}$$E25

$$m_{11}=m_1,m_{22}=m_2,m_{12}=c_1{m}_1+c_{2}m{}_2$$

$$m_i$$ is the ith particle’s mass. $$\overline{V_{ij}}$$ is the average molecular volume. Distribution function defines probability of finding particle at particular point $$r$$. In many literatures that have studied distribution function found it more versatile to use Laplace transform of this function $$G\left(s\right)$$.

$$G_{ij}\left(s\right)={\displaystyle {\int }_0^{\infty }r{g}_{ij}(r)e^{-sr}dr}$$E26

In this chapter the two formula which use RDF, we will encounter below integral equation that need expansion.

$$I={\displaystyle {\int }_{{\sigma }_0}^{\infty }r\varphi \left(r\right)rg\left(r\right)}dr={\sigma }_0^3{\displaystyle {\int }_1^{\infty }x\varphi \left(x\right)xg\left(x\right)dx=}{\sigma }_0^3\left\{{\displaystyle {\int }_0^{\infty }{x}^2\varphi \left(x\right)g\left(x\right)dx-{\displaystyle {\int }_0^1{x}^2\varphi \left(x\right)g\left(x\right)dx}}\right\}$$E27

On the right side of above equation from the right in the first equation we approximate distribution function with its values at contact points. This choice has been resulted from the behavior of molecules of which their repulsive interactions dominate their attractive potential. However, for the second term ($$I^{\prime }$$) we will use change in integrals to employ Laplace transform of RDF instead of RDF directly.

$$xg\left(x\right)=\frac{1}{2\pi i}{\displaystyle {\int }_{\gamma -i\infty }^{\gamma -i\infty }G\left(s\right)}{e}^{-sx}ds$$E28

Substituting above equation in $$I^{\prime }$$ we have

$$I^{\prime }={\displaystyle {\int }_0^{\infty }\phi \left(s\right)G\left(s\right)ds}$$E29

Where $$\phi \left(s\right)$$ represents

$$\phi \left(s\right)=\frac{1}{2\pi i}{\displaystyle {\int }_{\gamma -i\infty }^{\gamma -i\infty }x\varphi \left(x\right)}{e}^{-sx}dx$$E30

That indicates inverse Laplace of $$x\varphi \left(x\right)$$. So it suffices to just define inverse Laplace of potential function multiplied by $$x$$.

Therefore, Using Laplace transform of RDF $$G\left(s\right)$$ [9] quantum correction term for DY potential turn out to be

$$F^Q=\frac{h^2{N}_{A}n\beta }{24\pi }{\displaystyle \sum _{i,j}\frac{c_i{c}_j{\epsilon }_{ij}{A}_{ij}\overline{V_{ij}}}{m_{ij}{\sigma }_{ij}^0}\left({\lambda }_{ij}^2{e}^{{\lambda }_{ij}}G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)-{\upsilon }_{ij}^2{e}^{{\upsilon }_{ij}}G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)\right)}$$E31

$$c_i$$ is the i particle’s concentration and $$n$$ represents number density. $$A_{ij}$$, $${\lambda }_{ij}$$ and $${\upsilon }_{ij}$$ are controlling parameters of double Yukawa(DY). $${\epsilon }_{ij}$$ is the attractive well depth of mutual interacting potential.

3. Framework

The derivation of the thermodynamic and structural properties of a fluid system becomes a rather difficult problem when one wants to deal with realistic intermolecular interactions. For that reason, since the mid-20th century, simplifying attempts to (approximately) solve this problem have been devised, among which the perturbation theories of liquids have played a prominent role [10]. In this instance, the key idea is to express the actual potential in terms of a reference potential (that in terms of Ross perturbation theory Helmholtz free energy is expressed as of the “unperturbed” system) plus a correction term. This in turn implies that the thermodynamic and structural properties of the real system may be expressed in terms of those of the reference system which, of course, should be known. In the case of two component fluids, a natural choice for the reference system is the hard-sphere fluid, even for this simple system the thermodynamic and structural properties are known only approximately. Let us now consider a system defined by a pair interaction potential $$u(r)$$. The usual perturbation expansion for the Helmholtz free energy, $$F$$, to first order in $$\beta =1/k_{\beta }T$$, with $$T$$ being the absolute temperature and $$k_B$$ being the Boltzmann constant, leads to $$F$$. Common starting point of many thermodynamic perturbation theories is an expansion of the Helmholtz free energy, the resulting first-order prediction for a fluid composed of particles helium and hydrogen is given via the following equation.

$$F=F^t+F^Q+F^{HB}+F^{id}$$E32

The terms respectively are perturbation, Quantum, hard convex body and ideal terms. Perturbation term due to long range attraction of potential is given by [10]

$$F^t=2\pi n{\displaystyle \sum _{i,j}{c}_i{c}_j{\displaystyle {\int }_{{\sigma }_{ij}^0}^{\infty }{u}_{ij}(r)g_{ij}^{HS}(r,\rho ,{\sigma }_{ij})4\pi r^2\overline{V_{ij}}dr}}$$E33

Via Laplace transform of RDF (${rg}_{ij}^{HS}$) in calculation of first order perturbation contribution due to long-ranged attraction for DY potential we can employ below equation:

$$F^t=kT{\displaystyle \sum _{i,j}{c}_i{c}_j{\epsilon }_{ij}{\sigma }_{ij}^0{A}_{ij}{\overline{V}}_{ij}\left(e^{{\lambda }_{ij}}G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)-e^{{\upsilon }_{ij}}G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)\right)-\delta F^t}$$E34

$$\overline{V_{ij}}$$ the average molecular volume defined as:

$$\overline{V_{ij}}=1+\frac{({n^{\prime }}_i-1)}{{\sigma }_{ii}{}^3}\text{}\text{}\left[\frac{3}{2}\left({\sigma }_{ii}{}^2+{\sigma }_{ij}{}^2\right)l_i-\frac{1}{2}{\left(l_i\right)}^3-\frac{3}{2}{\left[{\left({\sigma }_{ii}+{\sigma }_{jj}\right)}^2-l_i{}^2\right]}^{\frac{1}{2}}{\sin }^{-1}\left[\frac{l_i}{{\sigma }_{ii}+{\sigma }_{jj}}\right]\left({\sigma }_{ij}{}^2\right)\right]$$E35

Where $${n^{\prime }}_i$$ define the number of element in a molecule, $$l_i$$ is distance of centre to centre for each molecule. ${\delta F}^t$ corresponds to the interval of $$\left[{\sigma }_{ij},{\sigma }_{ij}^0\right]$$ which long range attractive range is not further applicable. Consequently, we prefer to use the contact value of hard sphere RDF at $$r={\sigma }_{ij}$$. By this approach we can express this term as:

$$\delta F^t\approx \frac{n}{2kT}{\displaystyle \sum _{i,j}{c}_i{c}_j{\displaystyle {\int }_{{\sigma }_{ij}}^{{\sigma }_{ij}^0}{u}_{ij}(r)g_{ij}^{HS}({\sigma }_{ij})4\pi r^2\overline{V_{ij}}dr}}$$E36

$$g_{ij}^{HS}({\sigma }_{ij})$$ is the contact value of radial distribution function. $${\sigma }_{ij}$$ stands for separation distance at contact between the centers of two interacting fluid particles, with species $$i$$ and $$j$$. Although via minimization of Helmholtz free energy we can achieve value for hard sphere diameter, we preferred to use its analytical form due to its practical approach [17]. Hard sphere diameter will be calculated by means of Barker-Henderson equation as a function of interacting potential and temperature. Using Gauss-Legendre qudrature integration method we are able to evaluate its values numerically.

$${\sigma }_{ij}={\displaystyle {\int }_0^{{\sigma }_{ij}^0}\left(1-\exp \left(-\beta u_{ij}^{DY}(r)\right)\right)dr}$$E37

$$F^{HB},$$ Helmholtz free energy of hard convex body is given by following equation:

$$F^{HB}=a_{mix}\left(F^{HS}+F^{nd}\right)$$E38

Non-sphericity parameter $$a_{mix}$$ for the scaling theory [11] is defined as

$$a_{mix}=\frac{1}{3\pi }\frac{{\displaystyle \sum _{i,j}{c}_i{c}_j{V}_{ij}^{ef}{\left(V_{ij}^{eff}\right)}^{\prime }{\left(V_{ij}^{eff}\right)}^{″}}}{{\displaystyle \sum _{i,j}{c}_i{c}_j{V}_{ij}^{eff}}},\begin{array}{cc}& V_{ij}^{eff}=\frac{\pi }{6}{\sigma }_{ij}{}^3\end{array}\overline{V_{ij}}$$E39

$${\left(V_{ij}^{eff}\right)}^{\prime }$$ and $${\left(V_{ij}^{eff}\right)}^{″}$$ are the first and second partial derivatives of $$V_{ij}^{eff}$$ with respect to $${\sigma }_{ii}$$ and $${\sigma }_{jj}$$. From Boublik, Mansoori, Carnahan, Starling, Leland (BMCSL) [12, 8] with correction term of Barrio [13] on EOS, the Helmholtz free energy, $$F^{HS}$$ for hard sphere term becomes:

$$\frac{F^{HS}}{KT}=\frac{{\eta }_3\left[{\xi }_1+\left(2-{\eta }_3\right){\xi }_2\right]}{1-{\eta }_3}+\frac{{\eta }_3{\xi }_3}{{\left(1-{\eta }_3\right)}^2}+\left({\xi }_3+2{\xi }_2-1\right)\ln \left(1-{\eta }_3\right),$$

$${\xi }_1=\frac{3{\eta }_1{\eta }_2}{{\eta }_0{\eta }_3},{\xi }_2=\frac{{\eta }_1{\eta }_2}{{\eta }_3{}^2}\left({\eta }_4{z}_1+{\eta }_0{z}_2\right),{\xi }_3=\frac{{\eta }_2{}^3}{{\eta }_0{\eta }_3{}^3}$$

$$\begin{array}{l}{z}_1=2c_1{c}_2{\sigma }_{11}{\sigma }_{22}\left(\frac{{\sigma }_{11}-{\sigma }_{22}}{{\sigma }_{11}+{\sigma }_{22}}\right)\\ z_2=c_1{c}_2{\sigma }_{11}{\sigma }_{22}{}^3\left({\sigma }_{11}{}^2-{\sigma }_{22}{}^2\right)\end{array}$$E40

The correction term due to nonadditivity of the hard sphere diameter is the first order perturbation correction [14]

$$F^{nd}=-kT\pi n{c}_1{c}_2\left({\sigma }_{11}+{\sigma }_{22}\right)\left({\sigma }_{11}+{\sigma }_{22}-2{\sigma }_{12}\right)g_{12}^{HS}({\sigma }_{12})$$E41

In Eq. (41), $$g_{12}^{HS}({\sigma }_{12})$$ refer to as hard sphere radial distribution function at $$r={\sigma }_{12}$$ contact point by inclusion of Barrio and Solana correction on equation of state of BMCSL. Undoubtedly, the availability of the analytical HS RDF obtained from the solution to the corresponding Percus–Yevick (PY) equation represented a major step toward the successful application of the perturbation theory of liquids to more realistic inter-particle potentials. However, the lack of thermodynamic consistency between the virial and compressibility routes to the equation of state present in the PY approximation (as well as in other integral equation theories) is a drawback that may question the results derived from its use within a perturbation treatment. Fortunately, for our purposes, another analytical approximation for the RDF of the HS fluid, which avoids the thermodynamic consistency problem, has been more recently derived [15, 16]. We used improved RDF that yields exact asymptotic expression for the thermodynamic properties. However, we have used improved version of RDF that yields exact asymptotic expression for the thermodynamic properties. This have been derived by inclusion of Barrio and Solana correction on EOS of BMCSL at $$r={\sigma }_{12}$$ [9]

$$\begin{array}{l}{g}_{ij}^{HS}({\sigma }_{ij})=g_{ij}^{BMCSL}({\sigma }_{ij})+g_{ij}^{BS}({\sigma }_{ij})\\ g_{12}^{BMCSL}({\sigma }_{12})=\frac{1}{1-{\eta }_3}+\frac{3{\eta }_2}{{\left(1-{\eta }_3\right)}^2}\frac{{\sigma }_{ii}{\sigma }_{jj}}{{\sigma }_{ii}+{\sigma }_{jj}}+\frac{2{\eta }_2{}^2}{{\left(1-{\eta }_3\right)}^3}{\left(\frac{{\sigma }_{ii}{\sigma }_{jj}}{{\sigma }_{ii}+{\sigma }_{jj}}\right)}^2,\\ g_{12}^{BS}({\sigma }_{12})=\frac{1-{\delta }_{ij}{c}_i}{2}\frac{{\eta }_1{\eta }_2}{{\left(1-{\eta }_3\right)}^2}{\left(\frac{{\sigma }_{ii}{\sigma }_{jj}}{{\sigma }_{ii}+{\sigma }_{jj}}\right)}^2\left({\sigma }_{11}-{\sigma }_{22}\right)\left({\delta }_{ij}+\left(1-{\delta }_{ij}\right)\frac{{\sigma }_{22}}{{\sigma }_{11}}\right)\end{array}$$E42

$${\eta }_i=\frac{\pi }{6}n{\displaystyle \sum _j{c}_j{\sigma }_{jj}{}^i}$$

$${\delta }_{ij}$$ is the Kronecker delta function. For additive mixtures $${\sigma }_{ij}$$ is arithmetic mean of hard-core diameters of each species. Otherwise, the system is said to be non-additive.

The ideal free energy with $$N$$ particle for the atomic and molecular components of fluid mixture are given by,

$$F_{}^{id}(n,T)=\frac{3}{2}\ln \left(\frac{h^2}{2\pi kT{m}_1^{c_1}{m}_2^{c_2}}\right)+\ln n+{\displaystyle \sum _i{c}_i\ln c_i-1}$$E43

Compressibility factor of ideal term is one and $$Z^{HB}$$ would be estimated with the following derivation of related Helmholtz free energy

$$Z^{HB}=n\frac{\partial }{\partial n}\frac{F^{HB}}{kT}$$E44

For the perturbation term due to long rage attraction of potential tail employing (44) we will have

$$Z^t=\frac{2\pi n}{kT}{\displaystyle \sum _{i,j}{c}_i{c}_j{\epsilon }_{ij}{\sigma }_{ij}^0{A}_{ij}{\overline{V}}_{ij}\left(e^{{\lambda }_{ij}}\left(G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)-n\frac{\partial }{\partial n}G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)\right)-e^{{\upsilon }_{ij}}\left(G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)-n\frac{\partial }{\partial n}G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)\right)\right)-\delta Z^t}$$

$$\delta Z^t\approx \frac{n}{2kT}{\displaystyle \sum _{i,j}{c}_i{c}_j\left(g_{ij}^{HS}({\sigma }_{ij})+n\frac{\partial }{\partial n}{g}_{ij}^{HS}({\sigma }_{ij})\right){\displaystyle {\int }_{{\sigma }_{ij}}^{{\sigma }_{ij}^0}{u}_{ij}^{DY}(r)4\pi r^2\overline{V_{ij}}dr}}$$

Numerical integration has been used for calculation of $$\delta Z^t$$ in the range of $$\left[{\sigma }_{ij},{\sigma }_{ij}^0\right]$$. Expressions for first order perturbation and quantum correction term of compressibility factor are achievable via applying (44) for the free energy part of the quantum correction term.

$$Z^Q=\frac{h^2{N}_{A}n{\beta }^2}{24\pi }{\displaystyle \sum _{i,j}\frac{c_i{c}_j{\epsilon }_{ij}{A}_{ij}\overline{V_{ij}}}{m_{ij}{\sigma }_{ij}^0}\left({\lambda }_{ij}^2{e}^{{\lambda }_{ij}}\left(G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)-n\frac{\partial }{\partial n}G\left(\frac{{}^{{\lambda }_{ij}}}{{\sigma }_{ij}^0}\right)\right)-{\upsilon }_{ij}^2{e}^{{\upsilon }_{ij}}\left(G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)-n\frac{\partial }{\partial n}G\left(\frac{{}^{{\upsilon }_{ij}}}{{\sigma }_{ij}^0}\right)\right)\right)}$$E45

Summation over compressibility factors gives the total pressure of mixture

$$P=nkT\left(1+Z^{HB}+Z^t+Z^Q\right)$$E46

Defining Gibbs free energy provides information at critical points of phase stability diagram. Concavity and convexity of Gibbs diagram indicates if mixture is in one phase or not,

$$G=F+\frac{N}{n}P$$E47

Furthermore, Gibbs excess free energy is an appropriate measure in the definition of phase stability. Negative values for this energy describe stable state. This is expressed as

$$G_{xs}=G-{\displaystyle \sum _i{c}_i{G}_i^0}-NkT{\displaystyle \sum _i{c}_i\ln c_i}$$E48

That $$G_i^0$$ represents the Gibbs free energy of pure fluid of species i. Concentration-concentration structure factor is defined as

$$S_{cc}(0)=NkT{\left(\frac{{\partial }^2}{\partial c^2}G\right)}^{-1}{}_{T,P,N}$$E49

Compairing this equation with $$S_{cc}^{id}$$ enable us to define degree of hetero-coordination. In a given composition if $$S_{cc}(0)<<S_{cc}^{id}$$ then unlike atoms tend to pair as nearest neighbors (hetero-coordination) and when $$S_{cc}(0)>>S_{cc}^{id}$$ then like atoms are preferred as a neighbor.

3.1. Potentials

It is convenient to consider interacting potential with short-range sharply repulsive and longer-range attractive tail and treat them within a combined potential. The most practical method for the repulsive term of potential is the hard-sphere model with the benefit of preventing particles overlap. Furthermore, attractive or repulsive tails may be included using a perturbation theory. It is incontrovertible to generalize this potential to multi-component mixtures. This behavior is conveyed in double Yukawa (DY) potential which provides accurate thermodynamic properties of fluid in low temperatures and high density [18, 19]. At first we define DY potential as its effects on pressure of $$He-H_2$$ mixture has been studied in this work, written as:

$$u_{ij}^{DY}(r)={\epsilon }_{ij}{A}_{ij}\frac{{\sigma }_{ij}^0}{r}\left[e^{{\lambda }_{ij}\left(1-r/{\sigma }_{ij}^0\right)}-e^{{\upsilon }_{ij}\left(1-r/{\sigma }_{ij}^0\right)}\right]$$E50

$$A_{ij}$$, $${\lambda }_{ij}$$, $${\upsilon }_{ij}$$ are controlling parameters. These parameters for $$He$$ and $$H_2$$ are listed in table 1 with their reference [20].

	$$He-He$$	He-H2	H2-H2
$${\sigma }^0$$	2.634	2.970	2.978
$$A$$	2.548	2.801	3.179
$$\epsilon /k_B$$	10.57	15.50	36.40
$$\upsilon$$	3.336	3.386	3.211
$$\lambda$$	12.204	10.954	9.083

Table 1.

Potential parameters for $$He$$, $$H_2$$ interactions for DY potential [20].

For the atomic and molecular fluids studies in this mixture, these particles interact via a exponential six (exp-6) or Double Yukawa (DY) potential energy function [20]. The fluids considered in this work are binary mixtures that their constituents are spherical particles of two species, $$i$$ and $$j$$, interacting via pair potential $$u_{ij}(r)$$.

$$u_{ij}^{\exp -6}(r)=\{\begin{array}{l}\infty r<{\sigma }_{\infty }\\ {\epsilon }_{i,j}\frac{{\alpha }_{ij}}{{\alpha }_{ij}-6}\left(6\frac{6}{{\alpha }_{ij}}\exp ({\alpha }_{ij}(1-\frac{r}{{\sigma }_{\min ,ij}}))-{(\frac{{\sigma }_{\min ,ij}}{r})}^6\right)r>{\sigma }_{\infty }\end{array}$$E51

So we consider two-component fluid interacting via Buckingham potential $$u_{ij}(r)$$ between molecules of types $$i$$ and $$j$$. This potential is more realistic than square-well or Yukawa type potential for hydrogen isotope’s mixture [21] at high temperatures. Because of same atomic structure of hydrogen and its isotopes, the three constant of potential are same for hydrogen isotopes. These constants obtained experimentally from molecular scattering [22].

$${\sigma }_{\min }$$ indicate the range of interaction and the parameter $$\alpha$$ regulates the stiffness of repulsion. For hydrogen and helium type atoms these parameters have been organized in Table 2. It is well known that the long range attractive part of exp-6 potential is similar to Lenard–Jones potential.

In view of the energy equation (32), one can readily obtain equation for total pressure and different contributions to pressure from standard derivation of respective Helmholtz free energy. By the exp-6 potential, we have computed the Helmholtz free energy. The ten-point Gausses quadrature has been used to calculate integrals in quantum correction and perturbation contribution. The calculated pressure for $$D_2+T_2$$ fluid mixture with equal mole fraction and at temperature of $$T={100}^{o}K$$ is showed logarithmically in figure 1. As it is clear from this figure, the effect of hard sphere term of pressure in given rang of temperature is significant and the range of pressure variation is wider than ideal part. As it is mentioned earlier the difference between isotopes is simply related to the neutron number in each nucleus and correction due to difference in mass which involves in non-additive correction term doesn’t affect sensibly. In addition, this figure shows the predicted equimolar surface of the deuterium and tritium mixture for quantum correction term. This part is the most significant contribution at low temperature and varying smoothly in higher temperature. At very high densities, perturbation term contribution increases sharply with reducing density. Also, terms, $$P^t$$, $$P^Q$$, $$P^{HS}$$ and $$P^{id}$$, tend to infinity as $$\rho \to \infty$$.

	$$He-He$$	He-H2	H2-H2
$$\alpha$$	13.10	12.7	11.1
$$\epsilon /k_B$$	10.80	15.50	36.40
$${\sigma }_{\min }$$	0.29673	0.337	0.343

Table 2.

Potential parameters for He, H₂ interactions for exp-6 potential [20]

Figure 1.

Different contributions of pressure as a function of reduced density for $$T={100}^{o}K$$ for fluid mixture of deuterium and tritium

4. Results

For helium-hydrogen mixtures different parts of pressure due to correction terms and ideal parts have been showed in figure 2 at $$T={100}^{o}K$$. Ideal pressure at reduced density of approximately zero, to about 0.25 rises drastically. However, afterward it soars gently up to 100M (pa). Pressure due to hard sphere is the most significant contribution except that it is less than perturbation part at value of 1.5 for reduced density. Effects of perturbation and Quantum correction are important in high densities. In low densities, these contributions are insignificant and may possibly be ignored. Non-additive part has been caused by dissimilarity of particles which surges steadily from the beginning.

Figure 2.

Different contribution of correction terms on pressure of helium-hydrogen mixture at T=100, che=0.5 vs. reduce density

Gibbs excess free energy which is a measure for indicating phase stability of matters has been depicted in figure 3. Stability is limited to the areas that Gibbs excess free energy tends to negative values. This figure explains that stability rages for helium-hydrogen mixture at room temperature is confide in the boundaries in which helium concentration is less than 0.1.

Figure 3.

Gibbs excess free energy for helium-hydrogen mixture

Table 3 presents a comparison between results of pressure from this work using DY potential in place of exp-6, Monte–Carlo simulations and additionally study of reference [23] Obviously, there are appreciable adaption among our investigation results and MC which proves validity of our calculations. As Table 3, exhibits in low temperatures DY potential have more consistent results in comparison with exp-6. However, values of pressure extracted using DY potential cannot adjust with simulation resembling exp-6. Moreover, at higher temperatures after $$T={1000}^{o}K$$, DY potential is not good choice for evaluating EOS of hydrogen and helium mixture. We clarify our deduction presenting comparison between effects of these two potentials over wide ranges of temperatures.

$$T(K)$$	$$c_{He}$$	$${\rho }^{*}$$	$$\eta$$	$$P_{MC}$$	$$P$$[19]	$$P$$[23]	$$P_{DY}$$[19]
300	0.25	1.101	0.433	2.3090	2.7039	1.9664	2.8678
300	0.5	1.101	0.400	1.8560	1.7001	1.5729	1.8402
300	0.75	1.101	0.367	1.4240	1.2816	1.3160	1.3887
1000	0.5	1.223	0.335	4.5100	4.4205	4.1094	4.9406
1000	0.75	1.223	0.307	3.7150	3.5190	3.5904	3.9328
4000	0.5	1.376	0.247	12.4300	12.0832	12.1014	14.154
4000	0.5	1.572	0.282	16.3300	16.4485	16.4720	19.859

Table 3.

Comparison of values of pressure results from our study [19], Monte-Carlo simulation [24] and Isam Ali’s study [23].

Figure 4.

Comparision of efect of DY and EXP-6 potential on pressure of mixture in che=0.5, T=300 vs. Reduced density

Providing evidence of gradual divergence of DY and exp-6 potentials, a comparative figure has been made in figure 4 for helium-hydrogen mixtures. This figure shows more steepening effects of DY on total pressure of this mixture. Both potentials engender increase in pressure, except that, Buckingham affects moderately on pressure increase. The exp-6’s more steady behavior makes it adjustable with previous studies and MC simulation.

In figures 5, 6, 7, 8 we tried to give information about effects of quantum correction term on total pressure of helium-hydrogen and deuterium-tritium mixtures at the high reduced density of 1.3. This correction term has been plotted in 3-dimensional diagram in figure 5. This term is approximately zero for temperatures higher than 200 (K). Figure 6 represents that for hydrogen rich mixture at low temperature due to quantum effects pressure rise is significant. For effectual discussion on the effects of this term we have described $$P^Q/P$$ in figure 7 for helium-hydrogen and in figure 8 for deuterium tritium mixture. For the third picture increase in pressure is similar to what have been elaborated for figures 5 and 6. For figure 8 this manner remains analogous to helium-hydrogen mixture and temperatures next to 100 (K). However, for temperature lower than this it would behave inversely. For this range any increase in tritium concentration bears decrease in pressure.

Figure 5.

Pressure of quantum correction term at $$\rho s=1.3$$ for helium-hydrogen mixture.

Figure 6.

Total pressure from 50 K at ρ_S=1.3 for helium-hydrogen mixture.

Figure 7.

Fraction of quantum perturbation term to total pressure for helium-hydrogen mixture.

Figure 8.

Fraction of quantum perturbation term to total pressure for Deuterium-Tritium mixture.

5. Conclusion

An Equation of state of hydrogen–helium mixture has been studied up to 90G (pa) pressure and temperature equal to 4000◦K. We have used perturbation theory as an adequate theory for describing EOS of fluid mixtures. As well, by using this theory we can add extra distributive terms as perturb part which makes it more applicable than other theories. Considering this advantage, we can spread it out with additional terms for investigation on other states of matter like plasma in the direction of compares with experimental data. Otherwise, using simulation methods, for evaluating our theoretical results. Such as ab initio simulations with the code VASP,[25] which combines classical molecular dynamics simulation for the ions with electrons, behave in quantum mechanical system by means of finite temperature density functional theory [26]. In this chapter, two potentials have been presented, which we have used them for hydrogen isotopes and helium, and their mixtures. By means of comparison with Monte Carlo simulation and results of refrence [14] in Table 3 we could prove that exp-6 potential is more beneficial than DY in wider ranges of variables, since its application in this theory shows more convergent results in comparison with MC simulation [24]. Also exp-6 potential is a good choice of potential since it allows us to elevate temperature and density [28]. But as hydrogen molecules dissociation occurs [28] for pressures more than 100G (pa), this effect must be accounted. Therefore, we have restricted ourselves to pressures below 100G (pa).

Furthermore, we have used Wertheim RDF which enables us to use this EOS for extended values of temperature. As well, we have compared different contributions of pressure to represent which one is more effective in different density and temperature regimes. By finding the most effective parts of pressure contributions in each ranges of independent variables (Temperature, reduced density, mole fraction), we can omit the less significant parts which are considered ignorable in value, to decrease unnecessary efforts. Likewise, we can speculate from Fig. 1 that in low temperature and high densities, long range perturbation term has the most significant effect in comparison with other parts. On the other hand, hard sphere part can be assumed as the most noticeable part in high temperature ranges. Moreover, comparison of DY and exp-6 potentials effects, on pressure of this mixture has been studied to express benefits of using exp-6 potential for higher temperatures and densities. Additionally, as it is obvious in high temperature and density difference between effects of two potentials are considerable for this equimolar mixture. This discriminating property makes exp-6 potential preferable.

Furthermore, this approach has been used to evaluate EOS of $$D_2+T_2$$ mixture. Also, we have used this method to compare different contribution parts of pressure. These comparisons indicate that in low temperature quantum effects are more important, however in high temperatures, hard sphere part is the most effective. The last two three dimensional diagrams reveals the importance of quantum term in comparison with total pressure. However, for temperatures below 100 (k) for deuterium-tritium mixture negative pressure express that in low tritium concentration, deuterium rich fluid tend to consolidate.

6. Applications

One of the topics which can count on a great deal of interest from both theoretical and experimental physics is research in fluid mixture properties. These interests, not only comprise in the wide abundance of mixtures in our everyday life and in our universe but also the surprising new phenomena which were detected in the laboratories responsible for this increased attention. Mixtures, in general, have a much richer phase diagram than their pure constituents and various effects can be observed only in multi-component systems.

These kinds of studies have allowed a more complete modeling of mixture and consequently a better prediction and a more accurate calculation of thermodynamic quantities of mixture, such as activity coefficient, partial molar volume, phase behavior, local composition in general and have promoted a deeper understanding of the microscopic structure of mixtures.

Furthermore, for astronomical applications it is known that most of giant gas planets are like Jupiter is consisting primarily of hydrogen and helium. Modeling the interior of such planets requires an accurate equation of state for hydrogen-helium mixtures at high pressure and temperature conditions similar to those in planetary interiors [29]. Thus, the characterization of such system by statistical perturbation calculations will help us to answer questions concerning the inner structure of planets, their origin and evolution [29, 30].

In addition, in perturbation consideration of plasma via chemical picture, perturbation corrections will be included by means of additional free energy correction terms. Therefore, in considering transition behavior of molecular fluid to fully ionized plasma these terms are suitable in studying the neutral interaction parts. Consequently this will help us in studying inertial confinement fusion [31] and considering plasma as a fluid mixture in tokomak [32].