Convergence of Coupling-Parameter Expansion-Based Solutions to Ornstein–Zernike Equation in Liquid State Theory

Abstract

The objective of this paper is to investigate the convergence of coupling-parameter expansion-based solutions to the Ornstein–Zernike equation in liquid state theory. The analytically solved Baxter’s adhesive hard sphere model is analyzed first by using coupling-parameter expansion. It was found that the expansion provides accurate approximations to solutions—including the liquid-vapor phase diagram—in most parts of the phase plane. However, it fails to converge in the region where the model has only complex solutions. Similar analysis and results are obtained using analytical solutions within the mean spherical approximation for the hardcore Yukawa potential. However, numerical results indicate that the expansion converges in all regions in this model. Next, the convergence of the expansion is analyzed for the Lennard-Jones potential by using an accurate density-dependent bridge function in the closure relation. Numerical results are presented which show convergence of correlation functions, compressibility versus density profiles, etc., in the single as well as two-phase regions. Computed liquid-vapor phase diagrams, using two independent schemes employing the converged profiles, compare excellently with simulation data. The results obtained for the generalized Lennard-Jones potential, with varying repulsive exponent, also compare well with the simulation data. Solution-spaces and the bifurcation of the solutions of the Ornstein–Zernike equation that are relevant to coupling-parameter expansion are also briefly discussed. All of these results taken together establish the coupling-parameter expansion as a practical tool for studying single component fluid phases modeled via general pair-potentials.

1. Introduction

Liquid state theory deals with correlation functions, such as the pair-distribution function, to describe the structural and thermodynamics properties of fluids [1]. The correlation functions also form the basis of the modern theory of phase transitions in fluid systems [2]. The Ornstein–Zernike equation (OZE), which relates the direct and total correlation functions, is known to have excellent predictive power when supplemented with appropriate closure relations. The latter, as the name suggests, relates the two correlation functions and the interaction pair-potential. Although there are no exact closure relations, very accurate relations available in the literature provide results which compare excellently with those from computer simulations for homogeneous fluids [3]. There are also approaches wherein empirically introduced parameters in closure relations are tuned so as to impose thermodynamic consistency [4].

For realistic interaction potentials and accurate closure relations, the OZE is a highly non-linear integral equation, which is numerically solved using Newton’s method [5]. This approach works quite effectively for the case of purely repulsive interaction potentials [6]. However, a number of investigations show that the situation is quite different below the liquid-vapor phase transition temperature, which arises whenever there is an attractive component in the potential. Even simple models show the occurrence of multiple solutions and ‘no-real-solution-region’ in the liquid-vapor coexistence region [7]. This feature is most clearly illustrated in Baxter’s analytical solution of the adhesive hard sphere (AHS) model [8], which shows multiple solutions even with the diverging pair distribution function. The model employs the Percus–Yevick (PY) closure for hard spheres with adhesion. The proper solution branch from multiple solutions is easily selected so as to satisfy the ideal gas limit and spatial convergence of the pair distribution function. However, the model shows the emergence of a density domain in the two-phase region where all solutions are complex. It also shows the absence of a true spinodal line on the vapor side of the phase plane. Numerical methods also show similar features for the Lennard-Jones (LJ) potential and several important closure relations [9]. However, by using a different numerical technique with the PY closure, unique real solutions in the stable and metastable regions and multiple real solutions in the unstable (spinodal) region are also obtained [10]. It is difficult to enumerate all solution branches with any numerical scheme. On the other hand, analytical solutions within the mean spherical approximation (MSA) for the hardcore Yukawa (HCY) potential show the occurrence of two real and two complex solutions in the entire phase-plane. More specifically, it was found that there is a domain in the two-phase region where all solutions are complex. The model shows proper spinodal-lines on the vapor and liquid sides [11], as in van der Waals model and other mean field theories. Fortunately, the ‘no-real-solution-region’ in all models occurs inside the liquid-vapor coexistence dome, and careful strategies can be designed to map out the boundaries of the dome from numerical solutions of OZE [12].

An alternate approach to obtain useful approximations to fluid correlation functions is the thermodynamic perturbation theory (TPT) [13]. Here, the basic idea is to use the known (or accurately calculable) solutions for a reference potential and to incorporate the effects of perturbing potential via correction terms. Perturbation theories based purely on the fluid free energy mostly use the hard sphere model as the reference system [14]. This is because of the availability of simulation data on hard sphere systems in very accurate parameter forms [13]. While first order correction terms are easily computed by using the pair distribution function of the hard sphere system, three and four particle correlation functions are required even for second order term. Thus, the conventional TPT is restricted to low orders. The possibility of using OZE in higher order perturbation theory avoids the requirement of complicated correlation functions; however, it requires the computation of higher order derivative terms of the correlation functions [15]. A systematic development of this approach, called coupling-parameter expansion (CPE), to arbitrary high orders has provided good results for model systems in comparison to simulation data [16]. One of the most appropriate divisions of molecular potentials to reference and perturbing parts is based on the Weeks–Chandler–Andersen (WCA) prescription [17]. Recently, the CPE that employs this division and also incorporating density-dependent bridge functions is developed as a general method, making use of efficient Newton–Armijo and Krylov space-based solvers [18]. It has also been applied to compute phase diagrams of metallic fluids using effective pair-potentials derived from cohesive energy curves. However, there remains an important issue related to convergence of the CPE.

The main objective of the paper is to obtain results on the convergence of CPE, in particular, the two-phase region where OZE shows the occurrence of multiple and complex solutions. In what follows, CPE is first described briefly. Subsequently, its convergence is analyzed by using the analytical solutions for AHS and MSA-HCY systems. Important conclusions on convergence are obtained in this manner. Then, numerical results are obtained by using the general scheme for the LJ potential, which show convergence near the critical isotherm, and even in the two-phase region. A phase diagram of the LJ system is obtained via two independent schemes. Finally, the method is applied to generate the phase-diagrams for the general LJ(n,6) potentials and compared with simulation data. Important issues related to solution-spaces and bifurcation of solutions of OZE, which are relevant to coupling-parameter expansion, are also briefly discussed. The overall conclusion emerging from this analysis is that CPE is a practical approach for obtaining thermodynamic properties of one-component fluids interacting via general potentials. Thermodynamic models of fluid properties, particularly in the two-phase region, are essential in applications involving equations of state of liquid metals [19,20], which are important in high pressure physics. It is expected that CPE will find good use in generating the equations of state data [21] in the fluid domain.

2. Coupling-Parameter Expansion—General

The first step in introducing the general CPE [18] is to divide the inter-particle potential $U\left(r\right)$ into a reference (purely repulsive) and perturbation (purely attractive) components employing the WCA prescription [17]. These components are given by the following.

$$\begin{array}{ccc}\hfill U_0\left(r\right)& =& U\left(r\right)-U_{min},r\le r_{min},\hfill \\ & =& 0,r>r_{min},\hfill \\ \hfill U_a\left(r\right)& =& U_{min},r\le r_{min},\hfill \\ & =& U\left(r\right),r>r_{min}.\hfill \end{array}$$

(1)

Here, $U_{min}$ denotes the potential-minimum at position $r_{min}$. The purely repulsive part $U_0$ determines the structure of dense liquids while the attractive perturbing component $U_a$ gives rise to cohesion and phase changes [17]. When the coupling-parameter $\lambda$ is introduced into the potential, it reads as $U(r,\lambda )=U_0\left(r\right)+\lambda U_a\left(r\right)$. The magnitude of $\lambda$, which varies over $0\le \lambda \le 1$, determines the strength of the perturbation; thus, $U(r,0)$ is the reference part and $U(r,1)$ the full potential. Now, all correlation functions which determine the structure and thermodynamic properties of the fluid are also to be considered as functions of $\lambda$. In TPT, these functions are expanded as Taylor’s series in $\lambda$ about their reference part. For example, the pair distribution function $g(r,\lambda )$ is expanded as follows.

$$g(r,\lambda )=\sum _{n=0}^{\infty }\frac{1}{n!}{\lambda }^n{g}_n\left(r\right).$$

(2)

Here, $g_0\left(r\right)$ is the reference system’s pair distribution function and $g_n\left(r\right)$$(1\le n<\infty )$ denotes the $n^{th}$ derivatives of $g(r,\lambda )$ at $\lambda =0$. Similarly, Taylor’s series expansions are also assumed for other correlation functions: the direct correlation function $c(r,\lambda )$, the total correlation function $h(r,\lambda )=g(r,\lambda )-1$, and the indirect correlation function $y(r,\lambda )=h(r,\lambda )-c(r,\lambda )$. The subscript 0 in these functions denotes their corresponding functions for the reference system. The basic assumption in CPE is that these expansions converge and the series can be truncated at a suitable order. Thus, the main idea is that the contributions from the perturbing part of potential are accurately accounted for via the derivatives such as $g_n\left(r\right)$. While traditional perturbation theory is limited to the first two orders [13], the general formulation of CPE can incorporate quite higher order derivatives by taking recourse to the OZE.

The OZE defines the consistency relation between the short ranged direct correlation function $c(r,\lambda )$ and the long ranged total correlation function $h(r,\lambda )$. Similar to the pair-potential, both these functions tend to zero for large r. For one-component systems modeled using spherically symmetric potentials, the OZE is written as follows:

$$h(r,\lambda )=c(r,\lambda )+\rho \int c(|\overrightarrow{r}-\overrightarrow{r^{\prime }}|,\lambda )h(r^{\prime },\lambda )d{\overrightarrow{r}}^{\prime },$$

(3)

where $\rho$ is the number density of fluid particles. This is one equation in two unknowns, and it needs to be supplemented with a closure relation involving $h(r,\lambda )$ and $c(r,\lambda )$ and the pair-potential $U(r,\lambda )$. An exact closure is unfeasible; however, several approximate forms are available [1]. All the closures are generally expressed concisely in terms of a ‘bridge function’ $B(r,\lambda )$ in the following form:

$$c(r,\lambda )=exp[-\beta U(r,\lambda )+y(r,\lambda )+B(r,\lambda \left)\right]-y(r,\lambda )-1,$$

(4)

where $\beta ={\left(k_{B}T\right)}^{-1}$, $k_B$ is Boltzmann’s constant, and T denotes absolute temperature. Approximate forms of $B(r,\lambda )$ generate different closures; for instance, $B(r,\lambda )=ln[1+y(r,\lambda \left)\right]-y(r,\lambda )$ yields the PY closure while the hypernetted-chain closure (HNC) is obtained with $B(r,\lambda )=0$. A bridge function applicable to general potentials [22] is as follows:

$$\begin{array}{ccc}\hfill B(r,\lambda )& =& \sqrt{1+2\omega (r,\lambda )}-\omega (r,\lambda )-1,\hfill \\ \hfill \omega (r,\lambda )& =& y(r,\lambda )-\lambda {\rho }^{*}\beta U_a\left(r\right),\hfill \end{array}$$

(5)

which will be used throughout the present paper. Here and throughput the paper, ${\rho }^{*}=\rho {\sigma }^3$ is a dimensionless number density. This function has been shown to provide accurate thermodynamic properties for hard sphere as well as LJ systems.

2.1. Method of Solution

The solution of integral equation Equation (3) with a displacement kernel in Fourier space is easily expressed as follows:

$$\overline{y}(k,\lambda )=\frac{\rho {\overline{c}}^2(k,\lambda )}{1-\rho \overline{c}(k,\lambda )},$$

(6)

where the Fourier transform (denoted with a ‘bar’ throughout) and its inverse, for any function $q\left(r\right)$, are defined as the following.

$$\begin{array}{ccc}\hfill k\overline{q}\left(k\right)& =& 4\pi \underset{0}{\overset{\infty }{\int }}sin\left(kr\right)rq\left(r\right)dr,\hfill \\ \hfill rq\left(r\right)& =& {\left(2{\pi }^2\right)}^{-1}\underset{0}{\overset{\infty }{\int }}sin\left(kr\right)k\overline{q}\left(k\right)dk.\hfill \end{array}$$

(7)

The algebraic equations Equations (4) and (6), although the first is defined in r-space and the second in k-space, define a closed non-linear system within the integral equation formalism. Different implementations of Newton’s method to solve this system are now well developed [23]. However, they all have to deal with the problems mentioned earlier related to the ‘no-real-solution-region’ in the two-phase region.

In CPE, the Newton–Armijo solver is employed only for determining correlation functions of reference part. Starting with an initial guess $y_0$ (hence, $c_0$, ${\overline{c}}_0$, and ${\overline{y}}_0$), a correction $\Delta y$ is obtained by solving the following linear integral equation [18].

$$\Delta \overline{y}=\rho \frac{2{\overline{c}}_0+{\overline{y}}_0}{1-\rho {\overline{c}}_0}\overline{\left[\right(h_0+g_0{\Lambda }_0)\Delta y]}+[\rho \frac{{\overline{c}}_0^2}{1-\rho {\overline{c}}_0}-{\overline{y}}_0].$$

(8)

Here, ${\Lambda }_0={[1+2y_0]}^{-1/2}-1$ is the derivative $\partial B/\partial y$ of the bridge function $B(r,\lambda )$, defined in Equation (5), at $\lambda =0$. Fourier transforms and inverses are computed by using FFT (fast-Fourier-transform) algorithms. The linear integral equation implicitly defined in Equation (8) is best solved with Krylov space-based methods such as CGNR (Conjugate Gradient to Normal equations to minimize Residual) or GMRES (Generalized Minimum Residual). This is because it is possible to perform matrix-vector products via two Fourier transform operations. The Armijo scheme is used to limit the correction $\Delta \overline{y}$ to $\xi \Delta \overline{y}$ with $0<\xi <1$ in case error reduction does not occur at any iteration.

The derivatives of correlation functions are obtained from the solutions of the system of lineal integral equations.

$${\overline{y}}_n=[\overline{s_0}{}^2-1]\overline{\left[(h_0+g_0{\Lambda }_0)y_n\right]}+[\overline{s_0}{}^2-1]{\overline{w}}_n+{\overline{q}}_n,n\ge 1.$$

(9)

Here, ${\overline{s}}_0={[1-\rho {\overline{c}}_0]}^{-1}$ denotes the structure factor of the reference system. By denoting the derivatives of the bridge function also as $B_n={\Lambda }_0{y}_n+B_n^{*}$, the source terms, viz., $w_n$ and ${\overline{q}}_n$, are given by the following.

$$\begin{array}{ccc}\hfill w_n& =& -g_0\beta U_a{\delta }_{n1}+g_0{B}_n^{*}+\sum _{m=1}^{n-1}{C}_m^{n-1}[-\beta U_a{\delta }_{n-m,1}+y_{n-m}+B_{n-m}][y_m+c_m],n\ge 1,\hfill \\ \hfill {\overline{q}}_n& =& {\overline{s}}_0\sum _{m=1}^{n-1}{C}_m^n{\overline{c}}_{n-m}\rho [{\overline{y}}_m+{\overline{c}}_m],n\ge 1.\hfill \end{array}$$

(10)

Here, $C_m^n$ denotes the binomial coefficient. Note that the terms on right side involve only lower order derivatives, and, hence, are known at the n^th order. Furthermore, the linear operator in Equation (9) (for all orders) is identical to that in Equation (8). Finally, the derivatives $B_n^{*}$ are expressed as follows.

$$\begin{array}{ccc}\hfill B_n^{*}& =& -{\Lambda }_0{\rho }^{*}\beta U_a{\delta }_{n1}-(1+{\Lambda }_0)\sum _{m=0}^{n-2}{C}_m^{n-1}{\mathcal{O}}_{m+1}{\mathcal{O}}_{n-1-m},n\ge 1,\hfill \\ \hfill {\mathcal{O}}_n& =& (1+{\Lambda }_0)y_n+B_n^{*},n\ge 1.\hfill \end{array}$$

(11)

As the summation contributes only for $n\ge 2$, this equation provides $B_n^{*}$ explicitly at all orders. Derivations of these expressions are given in detail in the recent formulation of CPE mentioned earlier [18]. Concise algorithms for implementing the complete CPE method are also discussed there. (Note that the line $l=l$ in Algorithm A5 provided in this reference should be corrected as $l=2l$.)

2.2. Thermodynamic Functions

Once the derivatives are determined to the required order, explicit expressions for the correlation functions follow as functions of $\lambda$. Then, the series for $g(r,\lambda )$ provides an explicit expression for free energy F per particle [24].

$$\beta F(\rho ,T)=\beta F_0(\rho ,T)+\frac{\rho }{2}\int d\overrightarrow{r}\underset{0}{\overset{1}{\int }}d\lambda \beta U_a\left(r\right)g(r,\lambda ).$$

(12)

Here, $F_0$ is free energy per particle of the reference system. With the series expansion for $g(r,\lambda )$, the integral term over $\lambda$ is evaluated as ${\sum }_{n=0}^{\infty }{g}_n/(n+1)!$. By using the closure relation, which is also expressed as $g(r,\lambda )=exp[-\beta U(r,\lambda )+y(r,\lambda )+B(r,\lambda \left)\right]$, the above expression is recast as follows [24].

$$\begin{array}{ccc}\hfill \beta F(\rho ,T)& =& \beta F_0(\rho ,T)+\frac{\rho }{2}\int d\overrightarrow{r}[-g_a(r,1)+g(r,1)(y_a(r,1)+B_a(r,1))]\hfill \\ & -& \frac{\rho }{2}\int d\overrightarrow{r}\underset{0}{\overset{1}{\int }}d\lambda [y_a(r,\lambda )+B_a(r,\lambda )][\partial g_a(r,\lambda )/\partial \lambda ].\hfill \end{array}$$

(13)

Here, the subscript in $g_a(r,\lambda )$ (and similarly other functions) denotes the excess contribution over that of the reference potential, i.e., $g(r,\lambda )\equiv g_0\left(r\right)+g_a(r,\lambda )$. This form does not involve the attractive potential explicitly. Again, the integral term over $\lambda$ is obtained as ${\sum }_{n=1}^{\infty }{g}_n/(n-1)!{\sum }_{m=1}^{\infty }(y_m+B_m)/[(n+m)m!]$. Once $F(\rho ,T)$ is determined, pressure $\beta P={\rho }^2{(\partial F/\partial \rho )}_T$ and chemical potential $\beta \mu =\beta F+\beta P/\rho$ are readily computed. It is to be noted that this approach does not require integration over density or temperature.

The thermodynamic properties are also be obtained via the compressibility route. The (reduced) inverse compressibility ${\chi }^{-1}$ is provided in terms of a direct correlation function as follows.

$${\chi }^{-1}=\beta {(\partial P/\partial \rho )}_T=1-4\pi \rho \underset{0}{\overset{\infty }{\int }}c(r,1)r^{2}dr.$$

(14)

As the function $c(r,1)$ is short ranged, ${\chi }^{-1}$ is the most appropriate response function for computing pressure on an isotherm via integration with respect to $\rho$, whenever possible. Tje chemical potential $\mu$ is also determined by integrating the thermodynamic equation $\rho {(\partial \mu /\partial \rho )}_T={(\partial P/\partial \rho )}_T$, as follows.

$$\beta \mu =ln\left(\rho \right)+3ln\left({\lambda }_{th}\right)+\underset{0}{\overset{\rho }{\int }}\frac{d\rho }{\rho }({\chi }^{-1}-1).$$

(15)

The first two terms, where ${\lambda }_{th}=h/\sqrt{\left(2\pi m{k}_{B}T\right)}$ is the thermal wavelength, are the ideal gas contributions. This approach is always useful for computing $\beta F_0=\beta {\mu }_0-\beta P_0/\rho$ for the reference system and is employed in this paper. There are also explicit expressions for virial-pressure and chemical potential [12], involving $g(r,1)$, $c(r,1)$, and $B(r,1)$.

Having outlined the general algorithm for CPE, numerical investigations of convergence of the expansion for specific potentials are readily carried out. Before that, this analysis is conducted using analytical solutions of two model problems in the next two sections.

3. Coupling Parameter Expansion—AHS Model

Baxter’s AHS model deals with a system of hard spheres (of diameter $\sigma$) with an infinitely narrow attractive square well potential (of width $\Delta$) representing surface adhesion [8]. The well depth is taken as $\beta \varphi =ln[\tau \Delta /(\sigma +\Delta \left)\right]$, where the parameter $\tau$ represents an effective temperature, being zero at zero temperature and large at high temperatures. In the limit $\Delta \to 0$, this particular form reduces Mayer’s function to the form $f\left(r\right)$ = −1 + $(\sigma /12\tau )\delta (r-\sigma )$, where −1 represents the hard sphere contribution. By using virial expansion, it is also shown that the total correlation function reduces to $h\left(r\right)$ = −1 + $(\Gamma \sigma /12\tau )\delta (r-\sigma )$ for the range $0\le r\le \sigma$, where $\Gamma$ is a parameter to be determined employing OZE and PY closure. Baxter showed that ${\Gamma }^{*}$ = $\eta$ (1−$\eta$) $\Gamma$ satisfies the following quadratic equation:

$$\frac{1}{12}{\left({\Gamma }^{*}\right)}^2-[\tau (1-\eta )+\eta ]{\Gamma }^{*}+\eta (1+\eta /2)=0,$$

(16)

where $\eta =\pi \rho {\sigma }^3/6$ is the volume fraction of particles. This equation shows that the solutions are real in the entire range $0<\eta <1$ when $\tau >{\tau }_c=(2-\sqrt{2})/6\approx 0.09763$. However, they turn complex for a range ${\eta }_1\left(\tau \right)<\eta <{\eta }_2\left(\tau \right)$ around ${\eta }_c=(3/\sqrt{2}-2)\approx 0.1213$ when $\tau <{\tau }_c$. These boundaries are determined from the condition of vanishing discriminant in Equation (16), viz., ${[\tau (1-\eta )+\eta ]}^2-(\eta /3)(1+\eta /2)=0$. Thus the model predicts a gas-liquid phase transition with critical parameters ${\tau }_c$ and ${\eta }_c$. The proper root ${\Gamma }^{*}$ out of the two solutions that satisfy the ideal gas limit (for small $\eta$) is the negative root given by the following.

$${\Gamma }^{*}=6[\tau (1-\eta )+\eta ]-6{[{(\tau (1-\eta )+\eta )}^2-(\eta /3)(1+\eta /2)]}^{1/2}.$$

(17)

The analytical expression for pressure [8] that is obtained by integrating compressibility for $\tau >{\tau }_c$, and analytically continued to $\tau <{\tau }_c$ is given by the following.

$$\beta \frac{P}{\rho }=\frac{1+\eta +{\eta }^2}{{(1-\eta )}^3}-\frac{{\Gamma }^{*}(2+\eta )}{2{(1-\eta )}^3}+\frac{{\Gamma }^{*3}}{36\eta {(1-\eta )}^3}.$$

(18)

Here, the first term is the compressibility-pressure for hard spheres within the PY closure. Other terms involving ${\Gamma }^{*}$ are the contributions of surface adhesion. Similar expressions for chemical potential is also derived [25], and these expressions facilitate the determination of phase diagram [25] by using thermodynamic conditions. It is nice to note that the ‘no-real-solution-region’, where ${\Gamma }^{*}$ and, hence, P are undefined, is contained inside the co-existence region (see below). What is equally important is that the model shows the emergence of multiple solutions in all regions and a ‘no-real-solution-region’ in the two-phase region.

The equations of the AHS model are also derived as the lowest order solution to Baxter’s factorized form of OZE [26]. This analysis employed $ϵ=\Delta /(\sigma +\Delta )$ as a perturbation parameter [27] and directly showed that Equation (16) is accurate to $O\left(ϵ\right)$ and the parameters $\Gamma$ and $\tau$ are related to square well parameters as the following:

$$\tau ={\left(12ϵ\right)}^{-1}exp\left(\beta u_0\right),\Gamma =12\underset{1-ϵ}{\overset{1}{\int }}g\left[(\sigma +\Delta )x\right]xdx,$$

(19)

where $u_0<0$ is the depth of the potential. Now, all results of the model can be discussed in the usual phase-plane as the parameter $\tau$ is explicitly related to T.

The CPE is introduced by replacing the potential depth $u_0$ with $\lambda u_0$. Then, the functions $\tau \left(\lambda \right)$ and ${\Gamma }^{*}\left(\lambda \right)$, which is directly related to $\tau \left(\lambda \right)$, are expanded in a power series as the following.

$$\tau \left(\lambda \right)=\sum _{n=0}^{\infty }{\lambda }^n{\tau }_n,{\Gamma }^{*}\left(\lambda \right)=\sum _{n=0}^{\infty }{\lambda }^n{\Gamma }_n^{*}.$$

(20)

The series for $\tau \left(\lambda \right)$ is uniformly convergent with expansion coefficients ${\tau }_n={(12ϵn!)}^{-1}{\left[\beta u_0\right]}^n$ for $n\ge 0$. Substitution from Equations (20) into the quadratic equation, Equation (16), yields the following expansion coefficients ${\Gamma }_n^{*}$:

$$\begin{array}{ccc}\hfill {\Gamma }_n^{*}& =& \frac{{\omega }_n}{[{\tau }_0(1-\eta )-(1/6){\Gamma }_0^{*}+\eta ]},n\ge 1,\hfill \\ \hfill {\omega }_n& =& \frac{1}{12}\sum _{m=1}^{n-1}{\Gamma }_m^{*}{\Gamma }_{n-m}^{*}-(1-\eta )\sum _{m=1}^n{\tau }_m{\Gamma }_{n-m}^{*},n\ge 1,\hfill \end{array}$$

(21)

where ${\Gamma }_0^{*}$ satisfies an equation similar to Equation (16) but with ${\tau }_0$ in place of $\tau$.

$$\frac{1}{12}{\left({\Gamma }_0^{*}\right)}^2-[{\tau }_0(1-\eta )+\eta ]{\Gamma }_0^{*}+\eta (1+\eta /2)=0.$$

(22)

This equation is similar to the OZE for reference potential while Equation (21) is analogous to the solution of the system of linear integral equations defined in Equation (9). The parameters ${\tau }_0$ and ${\Gamma }_0^{*}$ correspond to the AHS model for depth $u_0=0$ but finite and small $ϵ$. This is taken as the reference system in lieu of the hard sphere system, because in order to obtain the hard sphere limit ${\tau }_0\to \infty$, it will be necessary to take the limit $ϵ\to 0$. The coefficients ${\Gamma }_n^{*}$ obtained up to sufficient order provides ${\Gamma }^{*}\left(1\right)$, which yields all thermodynamic properties within CPE. The series expansion for ${\Gamma }^{*}\left(\lambda \right)$ is also derivable from Equation (17). However, by writing the second square bracket term in this expression as $\sqrt{1-x}$, where $x=(\eta /3)(1+\eta /2)/{(\tau (1-\eta )+\eta )}^2$, this equation shows that the series fails to converge for ${\eta }_1\left(\tau \right)<\eta <{\eta }_2\left(\tau \right)$ when $\tau <{\tau }_c$, that is, within the ‘no-real-solution-region’.

The convergence properties just discussed are also elaborated in

Table 1. Convergence of CPE coefficients in AHS model.

AHS Model		${\mathsf{\Gamma }}_{\mathit{n}}^{*}/{\mathsf{\Gamma }}_{\mathit{n}-1}^{*}$	${\mathit{T}}^{*}$ = 0.45
n	$\mathbf{\eta }$= 0.05	$\mathbf{\eta }$= 0.15	$\mathbf{\eta }$= 0.25
1	2.11588	1.90467	1.69443
6	0.769139	1.93929	3.16261
11	0.823916	1.00009	0.764407
16	0.881579	0.948853	0.716068
21	0.90604	0.959392	0.842613
26	0.918982	0.971881	0.878875
31	0.927508	0.981006	0.880897
36	0.933696	0.987649	0.884413
41	0.938408	0.992688	0.889093
46	0.942117	0.996646	0.892829
51	0.945112	0.999839	0.895698

, where the ratio of the expansion coefficients $|{\Gamma }_n^{*}/{\Gamma }_{n-1}^{*}|$ versus n are given at three volume fractions for the subcritical temperature $T^{*}=0.45$ (see below). The ‘no-real-solution-region’ boundaries for this temperature are ${\eta }_1=0.06453$ and ${\eta }_2=0.1912$. The volume fractions $\eta =0.05$ and $0.25$ lie outside this region while $\eta =0.15$ lies inside. First and third columns illustrate convergence of CPE series, whereas the second column indicates a non-convergent series for $\lambda =1$. Thus, the real values of ${\Gamma }^{*}\left(1\right)$ obtained in the ‘no-real-solution-region’ at any order of truncation do not correspond to convergent series. However, this is not a serious issue, as will be discussed below.

The variation of ${\Gamma }^{*}\left(1\right)$ versus $\eta$ is compared with the exact solution (see Equation (17)) in Figure 1A (upper part) for two values of reduced temperature $T^{*}=k_{B}T/u_0$. These results correspond to the perturbation (width) parameter $ϵ=0.1$, which yields the critical temperature $T_c^{*}=0.4663$$({\tau }_c=0.09761)$. Graph of ${\Gamma }^{*}\left(1\right)$ versus $\eta$ obtained using CPE of 20th order $for{T}^{*}=0.5$$(\tau =0.1128)$ (curve-1) and exact results (symbols) compare excellently, and this case corresponds to $T^{*}>T_c^{*}$. Similar comparisons are found even at the 10th order. In the case of $T^{*}=0.4$$(\tau =0.06840)$, the range ${\eta }_1<\eta <{\eta }_2$ is the ‘no-real-solution-region’. Again, ${\Gamma }^{*}\left(1\right)$ versus $\eta$ from CPE of 20th order (curve-2) and the exact results (symbols) are compared. There are deviations between the two sets of results in the neighborhood of ${\eta }_1$ and ${\eta }_2$, otherwise they agree quite well. The portion of graph for ${\eta }_1<\eta <{\eta }_2$ obtained via CPE does not correspond to a convergent series. A subcritical isotherm for $T^{*}=0.45$$(\tau =0.09031)$ obtained via the compressibility route in CPE (curve-1) and the exact results (symbols) are compared in the lower part. Excellent agreement is found once again between the two sets of results, except in the ‘no-real-solution-region’ (blank portion). Similar agreement is obtained for chemical potential as well. Next, the phase diagram of AHS model is computed by using the analytical expressions for P and $\mu$ in the thermodynamic conditions, viz., $P_L=P_G$ and ${\mu }_L={\mu }_G$, where the subscripts L and G, respectively, denote the liquid and gaseous branches. The exact results for the phase diagram (blue line and symbols) are compared in Figure 1B with those obtained by using CPE of 50th order (curve-2) and 20th order (curve-3). These agree with one another, except for when it is close to the critical point. Critical temperature $T_c^{*}$ obtained at these orders are $0.4629$ and $0.4525$, respectively. The ‘no-real-solution-region’, lying within the dashed lines (curve-1), is also shown. These comparisons show that CPE provides useful results because the ‘no-real-solution-region’ where the expansion fails to converge lies completely inside the co-existence dome.

4. Coupling Parameter Expansion—MSA Model

The system of hard spheres (of diameter $\sigma$) attached with the attractive Yukawa potential is analytically solved within the MSA closure [28]. The Yukawa part, defined over the range $r\ge \sigma$, is expressed as $(-J\sigma /x)exp[-\zeta (r/\sigma -1\left)\right]$, where $\zeta$ and J are range and depth parameters, respectively. Analytical results of the model are useful for analyzing the convergence properties of CPE. Explicit expressions for all thermodynamic quantities exist in terms of a parameter $\Phi$; for instance, the reduced compressibility is obtained as follows:

$${\chi }^{-1}={(1-\eta )}^{-4}(1+2\eta +\Phi 12\eta {\zeta }^{-1}{e}^{\zeta }[(1+2\eta -6\eta /\zeta )(\mathcal{D}-1-\mathcal{D}{e}^{-\zeta }(1+\zeta ))+\mathcal{D}3\eta \zeta e^{-\zeta }]{)}^2,$$

(23)

where $\mathcal{D}$ is an explicit function of $\zeta$ and volume fraction $\eta$. Similar expressions are also derived for pressure computed via the virial and energy routes [11]. The central parameter $\Phi$ of the model satisfies a fourth degree polynomial equation of the following:

$$6{\eta }^2{\Phi }^4-X{\Phi }^3+12\eta K{\Phi }^2-KY\Phi +K^2=0,$$

(24)

where $K=\beta J$ and X and Y are specified functions of $\zeta$ and $\eta$. Note that this model has more solutions compared to the AHS model. This equation has two real and two complex solutions in the entire range of $\eta$ for K (proportional to $T^{-1}$) lower than a critical value $K^{*}\left(\zeta \right)$. Out of the two real solutions, the proper one must approach the ideal gas limit for small $\eta$ and should also vanish as $K\to 0$. For the region $K>K^{*}\left(\zeta \right)$, there is a ‘no-real-solution-region’ for a range ${\eta }_1\left(\zeta \right)<\eta <{\eta }_2\left(\zeta \right)$, because all solutions turn complex. Thus, the results of the model, although there are more solutions, are similar to those of the AHS model. However, an important difference is that the spinodal lines, which now occur on both the gaseous and liquid branches, encompass the ‘no-real-solution-region’. Thus, the occurrence of the latter may not be of much consequence in determining the phase diagram.

The CPE is introduced in the model by replacing K with $\lambda K$. Then the reference solution $(for\lambda =0)$ is that of hard spheres within the PY closure. As the proper real solution must approach zero as $\lambda \to 0$, the series expansion for $\Phi \left(\lambda \right)$ must start at first order and is given by the following.

$$\Phi \left(\lambda \right)=\sum _{n=1}^{\infty }{\lambda }^n{\varphi }_n.$$

(25)

Recursion relations for the expansion coefficients ${\varphi }_n$ readily follow on substituting this series in Equation (24):

$${\varphi }_1=\frac{K}{Y},{\varphi }_2=\frac{12\eta K{\psi }_{2,2}-X{\psi }_{3,3}}{KY},{\varphi }_n=\frac{12\eta K{\psi }_{2,n}-X{\psi }_{3,n+1}+6{\eta }^2{\psi }_{4,n+1}}{KY},n\ge 3,$$

(26)

where the additional coefficients ${\psi }_{2,n}$, ${\psi }_{3,n}$, and ${\psi }_{4,n}$ are obtained as the following.

$${\psi }_{2,n}=\sum _{m=1}^{n-1}{\varphi }_m{\varphi }_{n-m},n\ge 2;{\psi }_{3,n}=\sum _{m=1}^{n-2}{\varphi }_m{\psi }_{2,n-m},n\ge 3;{\psi }_{4,n}=\sum _{m=1}^{n-3}{\varphi }_m{\psi }_{3,n-m},n\ge 4.$$

(27)

After truncating the expansion at a suitable order, $\Phi \left(1\right)$ is to be substituted in the analytical expressions [11] to obtain the results within CPE. Note that CPE generates real solutions for all values of $\eta$ and K.

The convergence properties of CPE in the MSA-HCY model is different from that of the AHS model. The ratios of the expansion coefficients $|{\varphi }_n/{\varphi }_{n-1}|$ versus n are given in

Table 2. Convergence of CPE coefficients in the MSA-HCY model.

MSA-HCY Model		${\mathit{\varphi }}_{\mathit{n}}/{\mathit{\varphi }}_{\mathit{n}-1}$	K = 1.2
n	$\mathbf{\eta }=\mathbf{0.1}$	$\mathbf{\eta }=\mathbf{0.2}$	$\mathbf{\eta }=\mathbf{0.3}$
2	0.203451	0.233279	0.187291
6	0.568111	0.612217	0.454637
11	0.648161	0.689816	0.499086
16	0.678529	0.719187	0.51365
21	0.694581	0.734819	0.520534
26	0.704522	0.744575	0.524387
31	0.711287	0.751261	0.526771
36	0.71619	0.756136	0.528354
41	0.719906	0.75985	0.529465
46	0.722821	0.762777	0.53029
51	0.725169	0.765143	0.530938

for three volume fractions for $K=1.2$, which is in the subcritical range (see below). Volume fraction ${\eta }_1=0.1$ and $0.3$ lie outside the ‘no-real-solution-region’, while $\eta =0.2$ lies inside (see upper portion in Figure 2A). However, the data in this table illustrates convergence of CPE for all three $\eta$. Similar results are obtained for $K=1.4$ as well. Thus, CPE converges even the ‘no-real solution-region’ in this model.

The series solution $\Phi \left(1\right)$ versus $\eta$ is compared with exact solutions, which are simply the roots of Equation (24), for two values of the parameter K in Figure 2A. These results are for the range parameter $\zeta =2$, which yields the critical value $K^{*}=1.1354$. The lower part (for $K=0.8$) shows $\Phi \left(1\right)$, obtained from CPE of the 20th order (curve-1) and exact results (curve-2). The latter is the real root of Equation (24) discussed earlier. The two sets of results agree excellently, however, note that this case is in the region of K with two real solutions. The series solution converges quite fast, in fact, low order expansions (say 10th order) also yield similar comparisons. In the upper part of the figure (for $K=1.2>K^{*}$), the range ${\eta }_1<\eta <{\eta }_2$ is the ‘no-real-solution-region’. Graphs show $\Phi \left(1\right)$, obtained from CPE of 20th order (curve-1), and real (curve-2) and imaginary (curve-3) parts of the exact root. Note that CPE does not yield the complex solutions, and there are also significant deviations outside the ‘no-real-solution-region’. However, data in Table 2 have demonstrated convergence of CPE in all regions; thus, the graph for $\Phi \left(1\right)$ in the upper part also corresponds to a convergent series. Exact results for (reduced) inverse compressibility ${\chi }^{-1}$ (curve-2), obtained using analytical expressions [11], are compared in Figure 2B with those obtained by using CPE of 20th order (curve-1). The inset figure, which is for $K=1.4$, shows that ${\chi }^{-1}$ within CPE is non-negative in the entire range of $\eta$. Note that exact values of ${\chi }^{-1}$ are complex in the ‘no-real-solution-region’.

The results presented above for the AHS and MSA-HCY models show several features: (i) The CPE converges to exact results quite well in regions of phase plane where there are real solutions; (ii) the expansion fails to converge in ‘no-real-solution region’ of the AHS model, however, it converges even in that region for the MSA-HCY model; (iii) CPE is useful to derive the phase diagram by using thermodynamic conditions if it totally contains (as shown in the AHS model) the ‘no-real-solution region’.

5. Coupling Parameter Expansion—LJ Potential

This section presents numerical results for the more realistic LJ potential: $U\left(r\right)=4ϵ[{(\sigma /r)}^{12}-{(\sigma /r)}^6]$, where $ϵ$ and $\sigma$ are, respectively, the depth and range parameters. The density-dependent bridge function given in Equation (5) defines the closure relation for OZE. This is important because the predictive power of OZE depends solely on the accuracy of the bridge function. In obtaining numerical results, CPE (outlined earlier) is truncated at the 7th order, i.e., expansion including terms up to $g_6$ are used throughout the analysis. Reduced temperature $T^{*}=k_{B}T/ϵ$ and density ${\rho }^{*}=\rho {\sigma }^3$ are the variables used, and the argument $\lambda =1$ is omitted from different expressions.

The CPE around $\lambda =0$ provides real solutions throughout the phase plane if the series expansion converges. In the case of AHS and MSA models, as discussed in previous sections, convergence is clearly evident above the gas-liquid phase transition temperature. However, for the more realistic model considered here, convergence can only be tested numerically at somewhat lower orders of expansions. The LJ model (and bridge function) is also analyzed by applying Newton’s method to the full potential [9], where provisions to search for complex solutions are also builtin. This analysis reveals features similar to those found in AHS model, i.e., occurrence of a ‘no-real-solution-line’ on the gaseous side of an isotherm below critical temperate. Compressiblility is found to be finite but complex along this line, but the constant-volume heat capacity diverges. Although there is no true spinodal line on the gaseous side, the same analysis has found one to the left of the co-existence line on the liquid side. The emergence of complex isothermal compressibility within the HNC closure is also presented as a ‘fold-bifurcation’ [29] on isotherms for sub-critical temperatures. In fact, this feature is observed with several other closures [30], including the one considered here. However, it appears that the HNC closure does not yield a true critical point because the isothermal compressibility is finite on both boundaries (on gaseous and liquid sides) of the ‘no-real-solution region’ [4,31]. However, the series expansion in CPE about the repulsive component of the potential will not pick up these complex solutions. Moreover, as shown below, CPE is useful for computing the phase diagram, either via the compressibility or the free energy routes. The latter does not require data corresponding to phase points in the spinodal region.

First of all, it is necessary to check the accuracy of modeling the reference system resulting from WCA prescription. To that end, the pair distribution function $g_0(r/\sigma )$ versus scaled distance $(r/\sigma )$ (curve-1), which corresponds to $T^{*}=0.7290$ and ${\rho }^{*}=0.8446$, is compared against molecular dynamics simulation data (symbols) [32] in Figure 3A. Similar comparisons are shown in Figure 3B for the case $T^{*}=1.095$ and ${\rho }^{*}=0.6257$. Excellent agreement found in peak heights and variations in the entire range of $(r/\sigma )$ demonstrates the accuracy of the bridge function for the repulsive WCA potential. The compressibility factor $Z=\beta P/\rho$ versus the volume fraction, $\eta =\pi {\rho }^{*}/6$ is compared against simulation data (symbols) [33] in the inset figure for two temperatures: $T^{*}=0.75$ (top curve) and $T^{*}=2.0$ (bottom curve). Here, the agreement is also quite good, although slight differences are noted at higher $\eta$ for the low temperature case. All these results, obtained via the compressibility route, are found to be superior than those from the virial route.

Next, the convergence of the correlation functions $c(r/\sigma )$ and $g(r/\sigma )$ is analyzed in Figure 4 for single phase region and Figure 5 for two phase region. The former figure corresponds to the phase point ${\rho }^{*}=0.3$ and $T^{*}=1.45$, which is above the critical temperature, while the latter one is for ${\rho }^{*}=0.3$ and $T^{*}=1.15$ lying within the the spinodal region (see below). In both figures, the left and right panels (A and B) correspond to $c(r/\sigma )$ and $g(r/\sigma )$, respectively. The functions $c(r/\sigma )$ (curve-c), the reference function $c_0(r/\sigma )$ (curve-0), and six successive derivative terms $c_n(r/\sigma )/n!$ (curve-1 to curve-6) are shown in both left panels. Similarly, the two right panels show $g(r/\sigma )$ (curve-g), the reference function $g_0(r/\sigma )$ (curve-0), and six successive derivative terms $g_n(r/\sigma )/n!$ (curve-1 to curve-6).

Convergence to well defined functions, shown as $c(r/\sigma )$ and $g(r/\sigma )$, although based on numerical results, is evident at both phase points. Moreover, it is the first order term that contributes to the main correction. Note that the reference function $c_0(r/\sigma )$ is negative at both phase points; however, the graphs of $c(r/\sigma )$ for $r>\sigma$ have positive parts, which arise due to attractive interactions. This positive component is more pronounced at the spinodal region (Figure 5A), and it leads to negative compressibility at that point (see below). The reference functions $g_0(r/\sigma )$ have only small first peaks at both phase points. However, the first and second peaks in $g(r/\sigma )$ are more pronounced in the spinodal region (Figure 5B). This is because of lower temperature $T^{*}$ at this point, which amplifies the effects of attractive interactions. These profiles of $c(r/\sigma )$ and $g(r/\sigma )$ in Figure 5 (for the sub-critical case) are similar to those corresponding to one of the real solutions in the spinodal region [10], although they are obtained with PY closure. This point is discussed in more detail in Section 5.3 below.

Graphs of inverse compressibility ${\chi }^{-1}$ versus reduced density ${\rho }^{*}$, obtained using Equation (14) at successive orders, are shown in Figure 6. Here, the left panel (A) is for $T^{*}=1.45$ (in one-phase region), while the right panel (B) corresponds to $T^{*}=1.15$ (in spinodal region). Note that the graphs depict partial sums of the series terms. Thus, ${\chi }^{-1}$ were obtained using $c_0(r/\sigma )$ (curve-0), $c_0(r/\sigma )+c_1(r/\sigma )$ (curve-1), $c_0(r/\sigma )+c_1(r/\sigma )+(1/2)c_2(r/\sigma )$ (curve-2), etc., and finally $c(r/\sigma )$ (curve-6) are shown in the panels. Similarly to the case of correlation functions, the main contribution of attractive potential comes up at first order corrections, although it over-corrects a bit. Both panels show convergence of the curves to well defined profiles of ${\chi }^{-1}$. It is positive in the entire range of ${\rho }^{*}$ for $T^{*}=1.45$, as expected. However, for $T^{*}=1.15$, there is a region of ${\rho }^{*}$ where ${\chi }^{-1}$ is negative. Numerical values of ${\chi }^{-1}$ for three values of ${\rho }^{*}$, obtained with $c(r/\sigma )$, are also given in

Table 3. Partial sums of LJ model inverse compressibility series in CPE.

Order	Inverse Compressibility		${\mathit{T}}^{*}=1.15$	Inverse Compressibility		${\mathit{T}}^{*}=1.45$
	${\mathbf{\rho }}^{*}=\mathbf{0.091}$	${\mathbf{\rho }}^{*}=\mathbf{0.361}$	${\mathbf{\rho }}^{*}=\mathbf{0.631}$	${\mathbf{\rho }}^{*}=\mathbf{0.091}$	${\mathbf{\rho }}^{*}=\mathbf{0.361}$	${\mathbf{\rho }}^{*}=\mathbf{0.631}$
1	1.47757	4.53543	13.6363	1.46285	4.35925	12.6832
2	0.434523	−1.13238	2.19019	0.632217	−0.091533	3.77993
3	0.396093	−0.532527	2.79679	0.606776	0.283586	4.16096
4	0.369497	−0.491079	2.85451	0.592575	0.304848	4.19282
5	0.352317	−0.459612	2.88001	0.585209	0.318156	4.20459
6	0.342099	−0.43545	2.89254	0.581667	0.326561	4.20943
7	0.335684	−0.416063	2.89917	0.579865	0.332105	4.21159

. Data in column 3, for the case of ${\rho }^{*}=0.361$ and $T^{*}=1.15$, are for a phase-point in the spinodal region. The sub-critical isotherms, obtained by integrating ${\chi }^{-1}$ over density (starting from ${\rho }^{*}=0$), will show the familiar van der Waals loop due to the negative-value region. Thus, the thermodynamic properties obtained via CPE, although much more accurate, are similar to those of any first order mean field theory.

There is an issue related to the occurrence of negative compressibility (at any phase point) because it implies, by continuity requirements, that $1-\rho \overline{c}\left(k_0\right)$ would vanish at some positive value $k_0$. This is so because $\overline{c}\left(k\right)$ should tend to zero for sufficiently large k. Then, the transforms $\overline{y}\left(k\right)$ (see Equation (6)) and consequently $\overline{h}\left(k\right)$ would diverge at $k_0$. Thus, negative compressibility at any phase point implies a singularity at a positive $k_0$ in the Fourier transform $\overline{h}\left(k\right)$, which signals long ranged oscillations in the correlation function $h\left(r\right)$ [9]. Truncated Taylor’s series employed in CPE will not yield this singular function. In fact, this function is the solution to OZE with the series approximation to $\overline{c}\left(k\right)$. There is no inconsistency here, because it is known that OZE admits multiple solutions in the spinodal regions. As argued further in Section 5.3, only the series solution for $\overline{h}\left(k\right)$ obtained via CPE is relevant.

5.1. Phase Diagram of LJ Potential

Before obtaining the phase diagram, it is important to check the accuracy of CPE results on an isotherm close to the critical point. This is conducted in Figure 7A, which shows ${\chi }^{-1}$ versus ${\rho }^{*}$ (curve-1) on isotherm for $T^{*}=1.33$. The results of independent computations (symbols) are also shown [12]. The latter results are for identical models but are based on Newton’s method applied to the full potential, together with very careful strategies to mark out the spinodal region. The reduced temperature $T^{*}=1.33$ is just above the critical temperature $T_c^{*}=1.326$ and was obtained via CPE, as discussed next. The agreement between the two sets of results is quite good, thereby demonstrating convergence of ${\chi }^{-1}$ profile in CPE even close to the critical point. The inset figure shows pressure $\beta P$ (curve-1) and chemical potential $\beta \mu$ (curve-2) obtained via the compressibility route, i.e., by integrating ${\chi }^{-1}$ profile (see Equation (15)) from zero to ${\rho }^{*}$. Nearly flat regions in these graphs over a range of ${\rho }^{*}$ show the closeness of the isotherm to the critical isotherm.

Next, the phase diagram (co-existence curves) of the LJ system obtained via CPE is compared against simulation data (symbols) [34]. Two schemes are used for this purpose: (i) thermodynamic conditions (i.e., $P_L=P_G$ and ${\mu }_L={\mu }_G$ discussed earlier) via free energy route (curve-1) and (ii) Maxwell’s construction via compressibility route (curve-2). Note that the first scheme (see Equation (13)) is not dependent on the results within the spinodal region. The second scheme uses results in the negative compressibility region for obtaining P and its volume-integral on the liquid side of the phase-plane. The good agreement between the two schemes show that ${\chi }^{-1}$ profiles in the spinodal region are accurate interpolations between those in the gas and liquid sides. There are only slight differences between the results of the two schemes, along the liquid branch near the transition point. The simulation data agree better with the results of the former. The critical point parameters are also obtained using two methods: (i) solving the defining equations, ${\partial }_{v}P=0and{\partial }_v^{2}P=0$; and (ii) extrapolating values of ${\chi }_{min}^{-1}$ (over density) to zero on different isotherms in the one phase region. These methods provide almost identical critical point parameters: ${\rho }_c{\sigma }^3=0.299\left[0.312\right],k_B{T}_c/ϵ=1.326\left[1.316\right]$ and $P_c{\sigma }^3/ϵ=0.102\left[0.127\right]$, which compare well with simulation data given in square brackets. The spinodal lines on the gas and liquid sides (curve-3) are also shown in the figure.

The correlation functions and thermodynamic properties in the spinodal region, although converge to well defined profiles in CPE, are nonphysical, as they generate negative compressibility and van der Waals’s loop in pressure. However, their utility in computing the phase diagram (via thermodynamic conditions or Maxwell’s construction) is now demonstrated by comparing the results with those from independent methods. In any case, the weighted averaging of contributions of gas and liquid phases is needed to obtain physically acceptable properties.

5.2. Phase Diagrams of Generalized LJ(n,m) Potentials

The LJ potential is generalized to LJ($n,m$) family of potentials, with arbitrary power law exponents n and m, in the following expression.

$$U\left(r\right)=\frac{mϵ}{n-m}\left(\frac{n}{m}{)}^{n/(n-m)}\right[(\frac{\sigma }{r}{)}^n-(\frac{\sigma }{r}{)}^m].$$

(28)

The pre-factor is chosen so that potential depth is $ϵ$ for all n and m. Varying the exponent n, while keeping $m=6$, generates a family of potentials, also called Mie (n,6) potentials. The effective potential-range is altered on varying n: A larger n yields arrow and smaller n broad potentials. These potentials have applications in modeling systems in soft matter physics and have been investigated via molecular dynamics simulations [35] and integral equation theory [12].

The phase diagrams, obtained with CPE, are compared with simulation data (symbols) [35] in Figure 8A for potentials LJ(7,6) (curve-1), LJ(9,6) (curve-2), LJ(12,6) (curve-3), LJ(15,6) (curve-4), and LJ(20,6) (curve-5). These results are obtained via the compressibility route mentioned earlier. Perhaps the results obtained via CPE offer better comparison with simulation data than that obtained using Newton’s method applied to full potential [12]. It is also possible in CPE to approach the critical point closer on both branches. The gaseous branches and critical points compare well for all cases. However, there are differences on the liquid side for $n=7$, $n=9$, and $n=20$ cases. In fact, CPE results lie above simulation data for softer cases (smaller n) where the trend is just opposite for harder (larger n) cases. Results agree quite well on both branches for intermediate values $n=12$ and $n=15$. A possible reason for this is analyzed in Figure 8B, which shows the reduced isotherms for LJ(7,6) (curve-1) obtained via the virial (solid line), compressibility (dashed line), and free energy (dash-dot line) routes for $T^{*}=1$. Similar sets of graphs for potentials LJ(12,6) (curve-3) at $T^{*}=1$ and LJ(20,6) (curve-5) at $T^{*}=0.75$ are also shown in the panel. These results indicate that the bridge function [22] provides much better thermodynamic consistency for the intermediate cases ($n=12$ and 15) while important deviations are noted for softer and harder potentials. This might explain the differences in the phase diagrams on the liquid branch for the latter cases.

5.3. Solution-Spaces and Bifurcations of Solutions

The closure relation (Equation (4)) is solvable (in principle) so as to express $c\left(r\right)$ in terms of $h\left(r\right)$ explicitly. Then, OZE reduces to the following mapping: $h=\mathcal{O}h$, where $\mathcal{O}$ is a nonlinear operator. The specific details of $\mathcal{O}$ will depend on the closure relation, bride function, and inter-particle potential. However, the important point is that the solutions to OZE are simply the fixed-points of the operator $\mathcal{O}$, acting on functions defined in a suitable metric space [36]. Stability of the fixed-points and bifurcation of solutions, resulting due to variation of parameters (such as $\rho$ and T) in $\mathcal{O}$, will crucially depend on the solution-space chosen for $h\left(r\right)$ in addition to the specific properties of $\mathcal{O}$. If OZE is solved simply to obtain correlation functions, it is adequate to chose a solution-space $S_1$, wherein functions tend to zero for sufficiently large values of r. However, it is necessary to impose the constraint $I=\int h\left(r\right)d\overrightarrow{r}<\infty$, for employing $h\left(r\right)$ to compute thermodynamic quantifies, as reduced compressibility (along an isotherm) is given by $\chi =1+\rho I$. The integral constraint, which ensures the existence of 3D Fourier transforms of functions (used for computing structure factor), restricts the solutions to a sub-space $S_2$. Furthermore, numerical methods chosen to solve for the fixed-points, which are always of the iterative type, impose certain constraints on the solution so that iterations converge. It is common to modify the solution-algorithm, such as Newton’s method, and to look for complex solutions in regions where the method shows signs of diverging iterates [9].

The most efficient algorithms for solving OZE compute the convolution integral by using 3D Fourier transforms [6]. Baxter’s alternate formulation of OZE [8], applicable to cases where $c\left(r\right)=0$ for $r>R_0$ (cutoff range), is defined entirely in $0\le r\le R_0$. However, this formulation imposes the integral constraint mentioned above. The vanishing property of $c\left(r\right)$ is exactly valid within PY closure when potential $U\left(r\right)=0$ for $r>R_0$, as in the AHS model. Both these approaches seek solutions in the sub-space $S_2$; algorithms used with them are also modified to detect bifurcations from real to complex solutions at some $\rho$ along a sub-critical isotherm [7,9]. The region between such values of $\rho$ in the gas and liquid sides on the isotherm is termed ‘no-solution-region’, as the original aim was only to find real solutions. More detailed bifurcation analysis that was mentioned earlier [29], using the pseudoarc-length continuation (PAL) algorithm, has provided the complete bifurcation curve (with HNC closure). The same algorithm (PAL), as well as the modified Newton’s method, has also yielded real solutions within this region at somewhat lower temperatures, although not continuously across the ‘no-solution-region’. A numerical treatment of the AHS model [37] has shown similar features, wherein two pairs of complex solutions and one real solution are detected.

The OZE given in Equation (3) may also be written in the radial coordinate ($0\le r<\infty$) and solved numerically in a sufficiently long range $0\le r\le R$, together with builtin asymptotic conditions (for $r>R$) [10]. These solutions would then belong to the less restricted space $S_1$. Bifurcations of a real solution to three real solutions, along a sub-critical isotherm in the spinodal region, are found with this approach using PY closure [10]. Of the three solutions, two correspond to long ranged $h\left(r\right)$, which do not satisfy the integral constraint, and one to short ranged $h\left(r\right)$. Furthermore, the short ranged $h\left(r\right)$ and corresponding $c\left(r\right)$ yield negative compressibility along the part of the isotherm in the spinodal region. It is not known at present if this approach can also provide the complex solutions found in the sub-space $S_2$.

The real series obtained via CPE for $y\left(r\right)$ and, hence, $h\left(r\right)$ (in the stable, metastable, and unstable regions), provides short ranged solutions because the linear integral operator in Equation (9) involves only the correlation functions of the reference system. Furthermore, the source terms in this equation, for all orders, result from functions obtained at lower orders. Since the method uses 3D Fourier transforms, CPE seeks solutions in the sub-space $S_2$. The singular function $\overline{h}\left(k\right)$, obtained by substituting the series for $\overline{c}\left(k\right)$ in Equation (6), would be a first approximation to one of the long ranged $h\left(r\right)$; however, it lies outside the CPE method. This conclusion is based on the assumption that real solutions would be found in the spinodal region for the LJ model and closure used here, similarly to the PY case [10]. Finally, note that the closure relation rewritten as $g\left(r\right)=g_0\left(r\right)exp[-\beta U_a\left(r\right)+y_a\left(r\right)+B_a\left(r\right)]$ (see Equation (13)), together with the series for $y_a\left(r\right)$, can be used to obtain slightly improved results for pair correlation function.

6. Summary

The main aim in this paper is to analyze the convergence of CPE for solving the OZE. After providing brief details of CPE, its convergence is evaluated first for the two analytically solvable AHS and MSA models. Both these models show the emergence of a ‘no-real-solution-region’ in the two-phase domain of the phase plane. The CPE, which is based on series expansions around the correlation functions of the (repulsive) reference potential, always generate real solutions. Hence, it will not provide complex solutions to the OZE, obtained from several studies, in the two-phase region. The analytical arguments and numerical results presented for the AHS model show that CPE fails to converge in the ‘no-real-solution-region’. In fact, convergence is slow in the neighborhoods outside this region. However, this does not affect computation of the phase diagram of the AHS model because it completely contains the ‘no-real-solution-region’. Numerical results presented also show that CPE converges even in the ‘no-real-solution-region’ for the MSA-HCY model. The convergence issue is evaluated in detail for the LJ system modeled using an accurate density-dependent bridge function and 7th order solutions. Spatial profiles of correlation functions, corresponding to phase points above as well as below the critical point, are found to approach well defined functions. Similarly, density-profiles of inverse compressibility (above, below, and close to critical point) also converge toward appropriate functions. The phase diagrams of LJ model obtained via two schemes, one using the free energy and the other using the compressibilty routes, agree well and also with simulation data. The two schemes are independent as the free energy route uses only results corresponding to phase points outside the spinodal region. It is argued that CPE provides proper real solutions to OZE in the stable, metastable, and spinodal region. The whole procedure is also evaluated with reference to simulation data on the generalized LJ(n,6) family of potentials. These potentials with variable ranges have important applications in soft condensed matter. The main conclusion of this paper is that CPE provides a practical approach to evaluate the thermodynamic proprieties of one of the component fluids described by using general pair-potentials.