Equilibrium properties of penetrable soft spheres

Abstract

Using the Weeks-Chandler-Andersen separation scheme, we explored equilibrium properties of fully penetrable soft spheres with an attractive tail. When the radial distribution function of the reference system is computed by means of the integral equation theory using the hyper-netted-chain closure, this separation scheme predicts pressure-density isotherms, vapour-liquid phase coexistence, and surface tension in good agreement with the results from Monte Carlo simulations even though the attractive portion of the interparticle potential has a significant effect on the radial distribution function. Despite its simplicity, the model potential we studied, with a certain parameter value, exhibits a non-monotonic temperature dependence of the liquid phase density over a wide range of pressure.

GRAPHICAL ABSTRACT

Keywords:

• Penetrable soft spheres
• thermodynamic perturbation theory
• Weeks-Chandler-Andersen separation scheme
• hyper-netted-chain closure
• surface tension

1. Introduction

Soft spheres are characterised by a repulsive pairwise interaction that increases relatively slowly with decreasing interparticle distance when compared to commonly used functions such as the Lennard-Jones potential. Soft spheres find their utility in studying colloidal suspension and polymer solutions [1]. In recent years, soft spheres are often employed in meso scale simulations such as the dissipative particle dynamics (DPD) [2].

These particles are often fully penetrable in that the interparticle potential remains finite even when multiple particles occupy exactly the same position in space. The simplest of such bounded interparticle potential models is a step function, for which the potential energy is a positive constant up to some cut-off interparticle distance and vanishes beyond this distance. Schmidt developed a density functional for this model fluid based on the fundamental measure theory and predicted bulk fluid properties accurately at low packing fractions [3]. Other purely repulsive potentials, including Guassian core model [4] and exponential core model [5] were also studied.

The penetrable square well potential, a bounded version of square well potential, is the simplest model that exhibits vapour-liquid coexistence [6–8]. To study interfacial phenomena by means of DPD, Liu et al. [9] introduced a penetrable potential function constructed from cubic spline functions.

While particle based simulation techniques, such as molecular dynamics and Monte Carlo (MC) simulation [10], remain the primary tools for exploring the physical properties and dynamical behaviour of soft spheres, the computational demands imposed by these methods justify an effort to explore the same using a more analytical and computationally less demanding approaches. It is also of great interest to explore how well the existing theories of liquids [11] can predict the thermodynamic behaviour of soft spheres.

The integral equation theory is one such approach, in which one determines the radial distribution function by solving the Ornstein-Zernike (OZ) relation [11]. Various equilibrium properties of the system follows from this function. Since the OZ relation involves the total and the direct correlation functions, both of which are unknown, it needs to be supplemented by an additional so-called closure relation involving these functions.

For soft spheres, the hyper-netted-chain (HNC) closure [11] is known to yield a satisfactory result. The approach is particularly efficient computationally when applied to homogeneous phases. Recently, Malescio et al. [12, 13] used HNC to successfully reproduce phase diagrams in qualitative agreement with simulation and to predict the onset of mechanical instability of high density liquid phases.

However, HNC is not without some serious limitations. If it is applied to a system capable of exhibiting the vapour-liquid phase coexistence, the solution often becomes increasingly difficult to find as one approaches unstable (or even metastable) states. While this may be reasonable on a physical ground, it frustrates an attempt to predict the phase diagram accurately. It also implies that physical properties of an inhomogeneous system cannot be studied by means of a simple gradient theory [14–19], which requires the free energy density of a homogeneous system at all density values between those at the phase coexistence.

In this work, we consider three versions of HNC based approach. In the first version, which we shall refer to as HNC0, the full potential, consisting of both repulsive and attractive parts, will be treated by means of HNC. This version is subject to the limitations just indicated. In the second version, the potential energy is separated into a purely repulsive and a purely attractive parts according to the Weeks-Chandler-Andersen (WCA) scheme [20]. Only the repulsive part is treated using HNC. When the contribution from the attractive part is treated by the first order perturbation theory, this leads to a scheme we refer to as WCA/P. By changing how we treat the attractive part, we arrive at the third version to be referred to as WCA/V. The predictions of these approaches will be compared against the results from MC simulations.

We introduce the model potential in Section 2 and describe our theoretical approach for computing thermodynamic properties in Section 3. Section 4 briefly describes the details of MC simulations. The results of our computations are presented in Section 5. We conclude the article with a few remarks in Section 6. Appendix provides a short derivation of the key formula for calculating the surface tension.

2. Model potential

Following the work of Ref. [9], we consider the interparticle potential energy of the form (1) $\varphi (r)=\epsilon [W_{c}w(r/R_c)-W_{d}w(r/R_d)],$(1) where (2) $w(x)=\{\begin{array}{ll}1-6x^2+6x^3& (x\le 1/2)\\ 2(1-x{)}^3& (1/2<x\le 1)\\ 0& (x>1).\end{array}$(2) consists of cubic spline functions and has a continuous second derivative. In this work, we set $\epsilon =18.75$, $W_c=2$, $R_c=0.8$, $R_d=1$, and consider three values of $W_d$: $W_d=0.9$, 0.95, and 1. The resulting potential energy is shown in Figure 1.

Figure 1. Model potential given by Equations (1(1) $\varphi (r)=\epsilon [W_{c}w(r/R_c)-W_{d}w(r/R_d)],$(1) ) and (2(2) $w(x)=\{\begin{array}{ll}1-6x^2+6x^3& (x\le 1/2)\\ 2(1-x{)}^3& (1/2<x\le 1)\\ 0& (x>1).\end{array}$(2) ). $ϵ=18.75$, $W_c=2$, $R_c=0.8$, and $R_d=1$.

Using the WCA separation scheme [20], we split the potential energy at $r_{\min }$, where $\varphi (r)$ takes its minimum value ${\varphi }_{\min }$, into the repulsive part: (3) ${\varphi }_{\mathrm{r}}(r)=\{\begin{array}{ll}\varphi (r)-{\varphi }_{\min }& (r\le r_{\min })\\ 0& (r>r_{\min })\end{array}$(3) specifying the interparticle potential in the reference system, and the attractive part: (4) ${\varphi }_{\mathrm{a}}(r)=\{\begin{array}{ll}{\varphi }_{\min }& (r\le r_{\min })\\ \varphi (r)& (r>r_{\min }).\end{array}$(4)

3. Integral equation theory

Let $\psi :=\beta F/N$, where $\beta =1/k_{\mathrm{B}}T$ is the inverse temperature, F is the Helmholtz free energy, and N is the number of particles in the system. Then, (5) $\psi (T,\rho )=\ln ({\mathrm{\Lambda }}^3\rho )-1+{\psi }_{\mathrm{exc}}(T,\rho ),$(5) where Λ is the thermal wavelength of the particle and ρ is the number density of particles. The first two terms on the right-hand-side, taken together, represent the ideal gas contribution and ${\psi }_{\mathrm{exc}}$ is the contribution to ψ arising from the interparticle potential.

Under the WCA separation of the potential energy, ${\psi }_{\mathrm{exc}}$ is further divided into ${\psi }_{\mathrm{r}}$ and ${\psi }_{\mathrm{a}}$ due, respectively, to ${\varphi }_{\mathrm{r}}$ and ${\varphi }_{\mathrm{a}}$: (6) $\psi (T,\rho )=\ln ({\mathrm{\Lambda }}^3\rho )-1+{\psi }_{\mathrm{r}}(T,\rho )+{\psi }_{\mathrm{a}}(T,\rho ).$(6) For brevity and notational clarity, we present our method only for cases in which the WCA separation is employed (WCA/P and WCA/V). The formulation without this scheme (HNC0) can be recovered as a special case in which ${\varphi }_{\mathrm{r}}(r)\equiv \varphi (r)$ and ${\varphi }_{\mathrm{a}}(r)\equiv 0$. Accordingly ${\psi }_{\mathrm{r}}(T,\rho )\equiv {\psi }_{\mathrm{exc}}(T,\rho )$ and ${\psi }_{\mathrm{a}}(T,\rho )\equiv 0$ in HNC0.

3.1. Hyper-netted-chain closure

To determine ${\psi }_{\mathrm{r}}$, we find the radial distribution function $g_{\mathrm{r}}(r)$ of the reference system by solving the OZ relation (7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) combined with the HNC closure: (8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) where $c_{\mathrm{r}}(r)$ and $h_{\mathrm{r}}(r)=g_{\mathrm{r}}(r)-1$ are direct and total correlation functions, respectively.

Once $g_{\mathrm{r}}(r)$ is found, ${\psi }_{\mathrm{r}}$ is obtained using a thermodynamic identity (9) ${\psi }_{\mathrm{r}}(T,\rho )=-\frac{\beta p_{\mathrm{r}}(T,\rho )}{\rho }+\beta {\mu }_{\mathrm{r}}(T,\rho ),$(9) where $p_{\mathrm{r}}$ is the contribution to the pressure due to the repulsive interaction ${\varphi }_{\mathrm{r}}$ and is evaluated by the virial equation of state as (10) $-\frac{\beta p_{\mathrm{r}}(T,\rho )}{\rho }=\frac{1}{6}\beta \rho \int r\frac{\mathrm{d}{\varphi }_{\mathrm{r}}}{\mathrm{d}r}{g}_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(10) Under HNC [11], (11) $\begin{array}{rl}\beta {\mu }_{\mathrm{r}}(T,\rho )& =\frac{1}{2}\rho \int h_{\mathrm{r}}(r)[h_{\mathrm{r}}(r)-c_{\mathrm{r}}(r)]\mathrm{d}\mathbf{r}\\ & -\rho \int c_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.\end{array}$(11)

3.2. Perturbation theory

According to the first order perturbation theory, (12) ${\psi }_{\mathrm{a}}(T,\rho )=\frac{1}{2}\beta \rho \int {\varphi }_{\mathrm{a}}(r)g_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(12) Once ${\psi }_{\mathrm{a}}(T,\rho )$ is known, we can evaluate the contribution $p_{\mathrm{a}}$ to the pressure arising from ${\varphi }_{\mathrm{a}}$ by (13) $\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}={\left(\frac{\mathrm{\partial }{\psi }_{\mathrm{a}}}{\mathrm{\partial }\rho }\right)}_T.$(13) Alternatively, $p_{\mathrm{a}}$ may be evaluated using the virial equation of state as (14) $\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}=-\frac{1}{6}\beta \int r\frac{\mathrm{d}{\varphi }_{\mathrm{a}}}{\mathrm{d}r}{g}_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(14) Upon integration with respect to ρ, we obtain (15) ${\psi }_{\mathrm{a}}(T,\rho )={\int }_0^{\rho }\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}\mathrm{d}\rho .$(15) In the low density limit, $g_{\mathrm{r}}={\mathrm{e}}^{-\beta {\varphi }_{\mathrm{r}}}$. Thus, (16) $\underset{\rho \to 0}{lim}\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}=-\frac{2\pi }{3}\beta \int \frac{\mathrm{d}{\varphi }_{\mathrm{a}}}{\mathrm{d}r}{\mathrm{e}}^{-\beta {\varphi }_{\mathrm{r}}}{r}^3\mathrm{d}r.$(16) This may be evaluated by a numerical quadrature and gives the value of the integrand in Equation (15(15) ${\psi }_{\mathrm{a}}(T,\rho )={\int }_0^{\rho }\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}\mathrm{d}\rho .$(15) ) at $\rho =0$.

Equations (12(12) ${\psi }_{\mathrm{a}}(T,\rho )=\frac{1}{2}\beta \rho \int {\varphi }_{\mathrm{a}}(r)g_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(12) ) and (15(15) ${\psi }_{\mathrm{a}}(T,\rho )={\int }_0^{\rho }\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}\mathrm{d}\rho .$(15) ) define the methods we refer to as WCA/P and WCA/V, respectively. In both cases, the contribution from ${\varphi }_{\mathrm{a}}$ to the chemical potential is given by (17) $\beta {\mu }_{\mathrm{a}}(T,\rho )={\psi }_{\mathrm{a}}(T,\rho )+\frac{\beta p_{\mathrm{a}}(T,\rho )}{\rho }.$(17)

3.3. Phase coexistence

Densities of coexisting phases, ${\rho }_{\mathrm{v}}$ and ${\rho }_{\mathrm{l}}$, are determined by the equality of the pressures (18) $p(T,{\rho }_{\mathrm{v}})=p(T,{\rho }_{\mathrm{l}})$(18) and that of the chemical potentials: (19) $\mu (T,{\rho }_{\mathrm{v}})=\mu (T,{\rho }_{\mathrm{l}}).$(19)

3.4. Inhomogeneous system

We resort to the gradient theory to describe a vapour-liquid interface. In this theory, the free energy $F[\rho ]$ of an inhomogeneous system is a functional of the density profile $\rho (\mathbf{r})$ and is given by [14–19] (20) $\beta F[\rho ]={\int }_V[f_0(T,\rho )+f_2(T,\rho ){\left|\mathrm{\nabla }\rho \right|}^2]\mathrm{d}\mathbf{r},$(20) where $f_0:=\rho \psi$ is the free energy density of a bulk homogeneous phase and [11, 18, 19] (21) $f_2(T,\rho ):=\frac{1}{12}\int c(r)r^2\mathrm{d}\mathbf{r}.$(21) We used the random phase approximation to compute $c(r)$: (22) $c(r)=c_{\mathrm{r}}(r)-\beta {\varphi }_{\mathrm{a}}(r).$(22) The surface tension γ of a vapour-liquid interface is given by [19] (23) $\beta \gamma (T)=2{\int }_{{\rho }_{\mathrm{v}}}^{{\rho }_{\mathrm{l}}}\sqrt{f_2(T,\rho )[f_0(T,\rho )+\beta p-\beta \mu \rho ]}\mathrm{d}\rho ,$(23) where p and μ refer to the coexisting bulk phases. Appendix provides a short derivation of Equation (23(23) $\beta \gamma (T)=2{\int }_{{\rho }_{\mathrm{v}}}^{{\rho }_{\mathrm{l}}}\sqrt{f_2(T,\rho )[f_0(T,\rho )+\beta p-\beta \mu \rho ]}\mathrm{d}\rho ,$(23) ).

3.5. Numerical procedure

Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) were solved iteratively to find $c_{\mathrm{r}}(r)$ and $h_{\mathrm{r}}(r)$ at evenly spaced 501 grid points on $0\le r\le R_{\max }=5$. We examined the effect of doubling $R_{\max }$ or adding 500 more grid points under a few conditions. Except for a few instances to be stated explicitly, no discernible effect was observed.

At each temperature and density, HNC0 yields $p_{\mathrm{r}}$ and ${\mu }_{\mathrm{r}}$ directly without requiring an integration with respect to ρ. Accordingly, we have considered only those density values around the phase coexistence. Since WCA/V requires integration with respect to ρ as indicated by Equation (15(15) ${\psi }_{\mathrm{a}}(T,\rho )={\int }_0^{\rho }\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}\mathrm{d}\rho .$(15) ), we solved Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) for well over 200 density values. WCA/P shares the same set of data with WCA/V.

This allowed us to find an initial estimate for the coexisting densities by plotting μ versus p curve using ρ as a parameter and finding the density values at which the curve intersects itself. The initial estimate was further refined using the Newton-Raphson method [21].

In this and other computations, it becomes necessary to estimate various thermodynamic quantities and their density derivatives at arbitrary densities. For this purpose, we represented ${\psi }_{\mathrm{exc}}$, $\beta p_{\mathrm{exc}}:=\beta (p_{\mathrm{r}}+p_{\mathrm{a}})$, $\beta {\mu }_{\mathrm{exc}}:=\beta ({\mu }_{\mathrm{r}}+{\mu }_{\mathrm{a}})$, and $f_2$ using a cubic spline interpolant. In constructing the interpolant, we used the fact that the first three of these quantities are zero at $\rho =0$. To evaluate $f_2(T,\rho =0)$, we used the expression for $g_{\mathrm{r}}(r)$ in the low density limit in Equation (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ). Under the random phase approximation, this gives (24) $\underset{\rho \to 0}{lim}c(r)={\mathrm{e}}^{-\beta {\varphi }_{\mathrm{r}}(r)}-1-\beta {\varphi }_{\mathrm{a}}(r).$(24) We also need to specify the density derivatives of the above listed quantities both at the lowest ($\rho =0$) and at the highest densities for which data is available. Following a recommendation in Ref. [22], the derivative was estimated by representing the four data points at either end of the density values by a cubic polynomial.

A similar procedure was used for ${\psi }_{\mathrm{a}}$ in WCA/P so that $\beta p_{\mathrm{a}}/{\rho }^2$ can be evaluated by differentiation as indicated by Equation (13(13) $\frac{\beta p_{\mathrm{a}}(T,\rho )}{{\rho }^2}={\left(\frac{\mathrm{\partial }{\psi }_{\mathrm{a}}}{\mathrm{\partial }\rho }\right)}_T.$(13) ). In constructing the interpolant as described above, however, the following expression is available based on the low density approximation of $g_{\mathrm{r}}(r)$: (25) ${{\left(\frac{\mathrm{\partial }{\psi }_{\mathrm{a}}}{\mathrm{\partial }\rho }\right)}_T|}_{\rho =0}=2\pi \beta \int {\varphi }_{\mathrm{a}}{\mathrm{e}}^{-\beta {\varphi }_{\mathrm{r}}}{r}^2\mathrm{d}r.$(25)

4. Monte Carlo simulation

To assess the accuracy of our theoretical methods, we performed a series of MC simulations in canonical and Gibbs ensembles [23]. The bulk properties were computed not only for the system of particles interacting with the full potential φ, but also for the reference system in which the interparticle potential is purely repulsive and is given by Equation (3(3) ${\varphi }_{\mathrm{r}}(r)=\{\begin{array}{ll}\varphi (r)-{\varphi }_{\min }& (r\le r_{\min })\\ 0& (r>r_{\min })\end{array}$(3) ).

The canonical ensemble simulation for bulk properties were performed in a cubic box of various sizes chosen to contain a reasonable number of particles. For $\rho <1$, $1\le \rho <10$, and $10\le \rho \le 18$, the side length of the cube was 20, 12, and 8, respectively. The surface tension was evaluated in canonical ensemble simulations involving 10000 particles in a rectangular box of height 40 and a square base with the side length 8. At each state point considered, we performed 10 independent simulations, each involving sampling over $4\times {10}^5$ MC cycles. The resulting error bars were smaller than the symbols used to represent the data, and will not be shown.

The Gibbs ensemble (GEMC) simulation involved two cubic boxes, one for the vapour phase and the other for the liquid. We adjusted the total number of particles so that the average volume for these boxes are approximately ${12}^3$ and $8^3$, respectively. Only one simulation, involving sampling over $4\times {10}^5$ MC cycles, was performed at each temperature value.

In all cases, the periodic boundary condition was applied in each direction and equilibration was carried out over at least $4\times {10}^4$ MC cycles.

5. Results

To compare the HNC predictions, $g_{\mathrm{hnc}}(r)$, for the radial distribution function against the results, $g_{\mathrm{mc}}(r)$, of MC simulations, we quantified their difference by means of (26) $\Vert \mathrm{\Delta }g\Vert :=\frac{1}{R_g}{\int }_0^{R_g}|g_{\mathrm{hnc}}(r)-g_{\mathrm{mc}}(r)|\mathrm{d}r,$(26) where we dropped the usual factor of $4\pi r^2$ to emphasise the small r region ($r\le 1$), which is expected to have a more direct impact on the accuracy of the perturbation theory. Figure 2 shows $\Vert \mathrm{\Delta }g\Vert$ for $W_d=1$ and $R_g=R_{\max }=5$. Since the unweighted average of $g(r)$ over the interval $0\le r\le R_g$ is of the order of unity, we see that the error is generally well within a few % for both actual and reference systems. For the latter, there is a sudden increase in $\Vert \mathrm{\Delta }g\Vert$ at $\rho =1$. This is due to the finite size effect in MC simulations and results from our choice of the system volume described in Section 4. Using a larger system size reduced the magnitude of the jump, but only with a significant increase in the computational effort. For the actual system, the error is seen to grow rapidly towards the centre of the figure. This occurs as the system approaches the spinodal region (at $k_{\mathrm{B}}T=1$) or binodal lines (at $k_{\mathrm{B}}T=2$) either from the vapour or the liquid sides and the iterative solution of Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) becomes increasingly difficult. The error $\Vert \mathrm{\Delta }g\Vert$ in this region, especially for the liquid phases, was slightly larger with $R_{\max }=10$ than with $R_{\max }=5$ (using $R_g=5$ for both cases), presumably because correlations with longer wavelengths are permitted when $R_{\max }$ is larger.

Figure 2. Density dependence of the difference in the radial distribution functions between HNC and MC for the system interacting with the reference (${\varphi }_{\mathrm{r}}$) and the full (φ) potentials. $W_d=1$.

Figure 3 shows the p versus ρ isotherms for the actual system. A good agreement is observed between the MC simulation on the one hand and predictions from HNC0, WCA/P, and WCA/V on the other. We note that MC results for the liquid phase fall between the predictions of WCA/P and WCA/V. Isotherms from HNC0 agrees more closely with WCA/P except for low density liquid phases. We also note that the HNC0 results show a discontinuous jump in pressure (at $\rho \approx 10.3$ for $k_{\mathrm{B}}T=1$ and at $\rho \approx 8.4$ for $k_{\mathrm{B}}T=2$). This occurs at higher densities if $R_{\max }=10$. A convergent solution of Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) became impossible to find at these densities if $R_{\max }=20$. We took this as an indication that a homogeneous phase at these densities, at least according to HNC0, is unstable.

Figure 3. Pressure versus density isotherms. $W_d=1$.

When computing $p_{\mathrm{r}}$ by Equation (10(10) $-\frac{\beta p_{\mathrm{r}}(T,\rho )}{\rho }=\frac{1}{6}\beta \rho \int r\frac{\mathrm{d}{\varphi }_{\mathrm{r}}}{\mathrm{d}r}{g}_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(10) ), only the values of $g_{\mathrm{r}}(r)$ in the range $0\le r\le 1$ are involved. In Figure 4, we compare $g_{\mathrm{hnc}}(r)$ against $g_{\mathrm{mc}}(r)$. While HNC is seen to be reasonably accurate in predicting $g(r)$ for both actual and reference systems, deviation from $g_{\mathrm{mc}}(r)$ is noticeable for $0\le r\le 1$. In fact, when we set $R_g=1$ in Equation (26(26) $\Vert \mathrm{\Delta }g\Vert :=\frac{1}{R_g}{\int }_0^{R_g}|g_{\mathrm{hnc}}(r)-g_{\mathrm{mc}}(r)|\mathrm{d}r,$(26) ), $\Vert \mathrm{\Delta }g\Vert$ of the actual system increased approximately by a factor of 2 for the vapour and 5 for the liquid phases. Nevertheless, $\Vert \mathrm{\Delta }g\Vert$ is still very small, and hence the good agreement we observed for p between MC and HNC0 is not surprising.

Figure 4. The radial distribution function from HNC and an MC simulation. To improve visibility, graphs for $\rho =14$ and 18 are shifted upwardly by 1 and 2, respectively. $k_{\mathrm{B}}T=1$ and $W_d=1$.

A change in $\Vert \mathrm{\Delta }g\Vert$ of a similar magnitude was observed for the reference system when we used $R_g=1$ in Equation (26(26) $\Vert \mathrm{\Delta }g\Vert :=\frac{1}{R_g}{\int }_0^{R_g}|g_{\mathrm{hnc}}(r)-g_{\mathrm{mc}}(r)|\mathrm{d}r,$(26) ). In the case of WCA/P and WCA/V, however, the contribution from the attractive potential must also be included. Figure 4 reveals a considerable difference in the radial distribution function between the actual and the reference systems. This situation is in stark difference from a system of Lennard-Jones particles, in which $g(r)$ of the reference system (defined by the WCA scheme) and the actual system are very similar, with nearly perfectly coincident r values for peaks and valleys of $g(r)$.

To quantify the accuracy of the first order thermodynamic perturbation method, we evaluated the second order term using the approximate expression [11, 24]: (27) ${\psi }_{\mathrm{a}}^{(2)}(T,\rho )\approx -\frac{1}{4}\beta \rho {\left(\frac{\mathrm{\partial }{p}_{\mathrm{ref}}}{\mathrm{\partial }\rho }\right)}_T^{-1}\int {[{\varphi }_{\mathrm{a}}(r)]}^2{g}_{\mathrm{r}}(r)\mathrm{d}\mathbf{r},$(27) where $p_{\mathrm{ref}}:=\rho k_{\mathrm{B}}T+p_{\mathrm{r}}$. In Figure 5, we compare ${\psi }_{\mathrm{a}}^{(2)}$ against ${\psi }_{\mathrm{a}}$ given by Equation (12(12) ${\psi }_{\mathrm{a}}(T,\rho )=\frac{1}{2}\beta \rho \int {\varphi }_{\mathrm{a}}(r)g_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(12) ). We observe that these two quantities are comparable at gas phase densities, but the importance of ${\psi }_{\mathrm{a}}^{(2)}$ decreases rapidly with increasing density. At $\rho =10$, for example, ${\psi }_{\mathrm{a}}^{(2)}/{\psi }_{\mathrm{a}}$ is approximately 0.032 at $k_{\mathrm{B}}T=1.5$ and 0.01 at $k_{\mathrm{B}}T=0.5$.

Figure 5. The relative importance in the thermodynamic perturbation theory of the second order term ${\psi }_{\mathrm{a}}^{(2)}$ in comparison to the first order term ${\psi }_{\mathrm{a}}$ given by Equations (27(27) ${\psi }_{\mathrm{a}}^{(2)}(T,\rho )\approx -\frac{1}{4}\beta \rho {\left(\frac{\mathrm{\partial }{p}_{\mathrm{ref}}}{\mathrm{\partial }\rho }\right)}_T^{-1}\int {[{\varphi }_{\mathrm{a}}(r)]}^2{g}_{\mathrm{r}}(r)\mathrm{d}\mathbf{r},$(27) ) and (12(12) ${\psi }_{\mathrm{a}}(T,\rho )=\frac{1}{2}\beta \rho \int {\varphi }_{\mathrm{a}}(r)g_{\mathrm{r}}(r)\mathrm{d}\mathbf{r}.$(12) ), respectively. $W_d=1$.

The phase diagrams from WCA/P and WCA/V are also in good agreement with the results from GEMC simulations as shown in Figures 6 and 7. As with any other mean-field theory applied to model potentials with a harsh repulsive core, we observe that WCA/P and WCA/V both over-predict the critical temperature. For the range of $W_d$ values we considered, WCA/P and WCA/V predict very similar vapour phase densities, though the former tends to perform slightly better. In the case of the liquid phase, WCA/P is reasonably accurate for $W_d=0.9$ and 0.95. For $W_d=1$, however, GEMC results tend to interpolate WCA/P and WCA/V over a wide range of temperature. HNC0 results are generally in agreement with WCA/P at low temperatures. For the liquid phase, they approach WCA/V results as temperature is increased. At higher temperatures, HNC0 failed to provide a convergent solution for ρ values around the phase coexistence. As a result, the binodal lines terminated prematurely as seen in Figure 6.

Figure 6. Phase diagram showing vapour-liquid phase coexistence.

Figure 7. Vapor phase density at saturation.

Interestingly, MC, HNC0, WCA/P, and WCA/V all predict that the liquid phase density for $W_d=0.95$ increases with increasing temperature up to $k_{\mathrm{B}}T\approx 0.5$. This behaviour is related to the sign of the coefficient of thermal expansion, $-{\rho }^{-1}(\mathrm{\partial }\rho /\mathrm{\partial }T{)}_p$, of the liquid phase at and near the phase coexistence. To see this, we note that $\rho =\rho (T,p)$, and hence (28) $\mathrm{d}\rho ={\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }T}\right)}_p\mathrm{d}T+{\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}\right)}_T\mathrm{d}p.$(28) The Claussius-Clapeylon relation holds along the binodal lines: (29) $\mathrm{d}p=\frac{s_{\mathrm{v}}-s_{\mathrm{l}}}{1/{\rho }_{\mathrm{v}}-1/{\rho }_{\mathrm{l}}}\mathrm{d}T,$(29) where $s_{\mathrm{l}}$ and $s_{\mathrm{v}}$ denote the entropy per particle in the liquid and the vapour phases, respectively. Eliminating dp from Equation (28(28) $\mathrm{d}\rho ={\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }T}\right)}_p\mathrm{d}T+{\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}\right)}_T\mathrm{d}p.$(28) ), we arrive at (30) ${\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }T}\right)}_{\mathrm{vle}}={\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }T}\right)}_p+\frac{s_{\mathrm{v}}-s_{\mathrm{l}}}{1/{\rho }_{\mathrm{v}}-1/{\rho }_{\mathrm{l}}}{\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}\right)}_T,$(30) where the subscript ‘vle’ indicates that the derivative is taken while maintaining the vapour-liquid phase coexistence. Insofar as the second term on the right hand side is expected to be positive, $(\mathrm{\partial }\rho /\mathrm{\partial }T{)}_{\mathrm{vle}}>0$ if $(\mathrm{\partial }\rho /\mathrm{\partial }T{)}_p>0$. Figure 8 shows that the latter derivative from MC simulations indeed is positive up to $k_{\mathrm{B}}T\approx 0.4$ over a very wide range of pressure. HNC0, WCA/P, and WCA/V all correctly capture this behaviour at least qualitatively.

Figure 8. Temperature dependence of the the liquid phase density at p = 0.01, p = 1, and p = 10. $W_d=0.95$. For each method and at each temperature, the density is larger for a higher pressure.

As shown in Figure 9, WCA/P and WCA/V both predict the surface tension γ reasonably well. WCA/V is more accurate at $W_d=1$ but WCA/P becomes comparably accurate as $W_d$ is decreased. For all $W_d$ values we studied, WCA/P and WCA/V over-predict γ at very high temperature values. This is a consequence of these theories over-predicting the critical temperature. We recall that HNC0 fails when applied to unstable phases and hence cannot be used to predict γ within the framework of the gradient theory.

Figure 9. Surface tension γ of the vapour-liquid interface at saturation.

For a sufficiently strong attractive potential, particles in the system may collapse to form a small blob. According to the Fisher and Ruelle criterion [25, 26], the system at $k_{\mathrm{B}}T=0$ is unstable with respect to this collapse if (31) $\int \varphi (r)\mathrm{d}\mathbf{r}<0,$(31) which gives $W_d>1.024$ for the model potential we are considering in this work. At non-zero temperatures, the entropic effect is expected to prevent the collapse if $W_d$ is exactly 1.024. As $W_d$ is increased, however, the collapse should be observable even at non-zero temperatures.

Taking the failure of the iterative solution of Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) as indicating the onset of this instability, Malescio et al. [12, 13] showed that HNC0 predicts the onset in an excellent agreement with the Fisher and Ruelle criterion for several soft sphere potentials.

In contrast, we do not expect WCA/P or WCA/V to be capable of predicting this form of instability since the reference system for which HNC equation is solved is purely repulsive. Nevertheless, we found that WCA/V (as well as HNC0) predicts a form of mechanical instability (in addition to the spinodal) when the density is increased far beyond the liquid density at the vapour-liquid coexistence.

As an example, we considered three $W_d$ values: $W_d=1.020$, 1.024, and 1.028. The onset of this mechanical instability is identified by (32) ${\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }\rho }\right)}_T=0.$(32) This equation is solved at various temperatures and the solution, ${\rho }_{\star }$, is shown in Figure 10. At least for $W_d=1.020$ and 1.024, WCA/V is seen to considerably underestimate ${\rho }_{\star }$ compared to HNC0. In both HNC0 and WCA/V, p decreased rapidly when ρ is increased beyond ${\rho }_{\star }$. The iterative solution of Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) eventually failed to converge at ${\rho }_{\max }$.

Figure 10. The density ${\rho }_{\star }$, beyond which a high density fluid phase is mechanically unstable, plotted versus temperature. An open symbol indicates a fluid phase with a negative pressure.

In Figure 10, open symbols indicate that $p(T,{\rho }_{\star })<0$. Since the pressure takes the local minimum at the liquid phase spinodal density ${\rho }_{\mathrm{l}}^{\mathrm{sp}}$ before it reaches the local maximum at ${\rho }_{\star }$, we have $p(T,{\rho }_{\mathrm{l}}^{\mathrm{sp}})<p(T,{\rho }_{\star })<0$. The vapour phase, if existed, should have a positive pressure. As a result, the system has no vapour-liquid coexistence at the temperature and the $W_d$ value indicated by an open symbol in Figure 10. This means that, as the temperature is dicreased, the binodal lines terminate at T satisfying $p(T,{\rho }_{\star })=0$, where ${\rho }_{\star }$ for a given $W_d$ is a function of T only.

This is illustrated in Figure 11, which shows the liquid densities at saturation, ${\rho }_{\star }$, and ${\rho }_{\max }$ for $W_d=1.024$. Compared to HNC0, WCA/V predicts the termination of the binodal line (due to $p(T,{\rho }_{\star })$ being zero) at a much higher temperature. The HNC0 predictions follow the GEMC results more closely in this temperature range.

Figure 11. Liquid phase density at saturation (${\rho }_{\mathrm{l}}^{\mathrm{eq}}$) shown with the onset of mechanical instability (${\rho }_{\star }$) and the maximum density (${\rho }_{\max }$) above which iterative solution of Equations (7(7) $h_{\mathrm{r}}(r)=c_{\mathrm{r}}(r)+\rho \int c_{\mathrm{r}}(|\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }}|)h_{\mathrm{r}}(r^{\mathrm{\prime }})\mathrm{d}{\mathbf{r}}^{\mathrm{\prime }}$(7) ) and (8(8) $c_{\mathrm{r}}(r)=h_{\mathrm{r}}(r)-\ln [h_{\mathrm{r}}(r)+1]-\beta {\varphi }_{\mathrm{r}}(r),$(8) ) fails. $W_d=1.024$. An open symbol for ${\rho }_{\star }$ indicates a negative pressure.

However, the high density phase in GEMC simulation crystallised at and below $k_{\mathrm{B}}T=0.78$, thus preventing a meaningful comparison between HNC0 and GEMC. The crystal phase had 4- and 6-fold rotational symmetries, pointing to the face-centred-cubic structure. As with ${\rho }_{\star }$, WCA/V predictions of ${\rho }_{\max }$ are considerably smaller than those of HNC0. In agreement with what we saw for $W_d=0.9$, 0.95, and 1, the binodal line from HNC0 terminated prematurely as the temperature was increased.

WCA/P exhibited a rather peculiar behaviour at $k_{\mathrm{B}}T=0.2$ and 0.4. The p versus ρ isotherms have a local maximum at ${\rho }_{\star }(>{\rho }_{\mathrm{l}}^{\mathrm{sp}})$, which is followed by a local minimum. With a further increase in ρ, p increased rapidly before the iterative procedure failed at the same ${\rho }_{\max }$ value as for WCA/V. (WCA/P and WCA/V share the same value for ${\rho }_{\max }$, but not for ${\rho }_{\star }$). The two extrema eventually merged into a single inflection point, which disappeared at around $k_{\mathrm{B}}T=0.6$. The sharp increase in p and the eventual failure of the iterative process were observed at all temperature values.

A failure of the iterative process is often identified as an onset of some instability. To gain an insight into the nature of this instability, we examined the radial distribution function and noticed an emergence and a sudden growth of a peak near r = 0. This is shown in Figure 12 for $W_d=1.024$ using the value of $g_{\mathrm{hnc}}(r)$ at r = 0.01, this value of r being the smallest available from our computation. Similar plots for other values of $W_d$ (1.020 and 1.028) are omitted in the figures as they are indistinguishable from the graphs for $W_d=1.024$. The figure clearly indicates a significant overlapping of particles with increasing ρ as envisaged in the literature [12, 13, 26].

Figure 12. Density dependence of the radial distribution function at r = 0.01 using the HNC closure. $k_{\mathrm{B}}T=1$ and $W_d=1.024$. The filled circules indicate the values at the onset of mechanical instability ($\rho ={\rho }_{\star }$).

The rapid increase in $g_{\mathrm{hnc}}(r=0.01)$ occurs at lower densities for the reference system than for the actual system. This is expected and leads to the under-prediction of ${\rho }_{\star }$ and ${\rho }_{\max }$ by WCA/V mentioned earlier: The energetically favourable interaction is operational only for $r_0\le r\le R_d=1$, where $r_0$ is determined by $\varphi (r_0)=0$, and this delays the overlapping of the particles in the actual system.

We emphasise that the overlapping of particles is observed even in the reference system, in which the interparticle potential is non-negative. There appears to be no reason for such particles to condense into a small blob. Even though a repulsive force from a single particle may be weak, a strong repulsive force would be generated if multiple of such particles fully overlap. Thus, it does not appear unreasonable to speculate that, at least in the case of the reference system and by extension in the actual system with small $W_d$ values, the new phase may be a crystal, in which each lattice point consists of multiple overlapping particles. A systematic exploration of this possibility would require a density functional theory of crystals [27–29] along with extensive and very demanding (due to the very high density values involved) MC simulations, which is beyond the scope of this study.

6. Concluding remarks

We examined the accuracy of the WCA separation scheme for fully penetrable soft spheres with a bounded potential energy. When the total and direct correlation functions of the reference system are determined by the HNC closure, this separation scheme leads to predictions of equilibrium properties of both homogeneous and inhomogeneous systems in a good agreement with MC simulations. In contrast to more commonly studied potentials characterised by a harsh repulsive core, however, the radial distribution function of the reference system differs noticeably from that of the actual system.

The model potential we studied, despite its simplicity, leads to a non-trivial equilibrium behaviour. For example, the liquid phase density, as a function of temperature, is not monotonic over a very wide range of pressure at least for one of the $W_d$ values we examined. With increasing density, the liquid phase eventually becomes mechanically unstable. If the attractive part of the potential energy is sufficiently strong, we expect the liquid phase to collapse and form a high density blob. However, the reference system (without an attractive potential) is likely to form a crystal phase instead, in which each lattice point consists of multiple overlapping particles. Thus, unless the attractive potential is sufficiently strong, the above mentioned mechanical instability might point to a crystal phase formation. A systematic exploration of this possibility, however, is left to a future work.

Acknowledgments

Computations reported here were made possible by a resource grant from the Ohio Supercomputer Center [30]. We are grateful for helpful comments from Professor Lisa M. Hall on the integral equation theory during the early stage of this work. A preliminary GEMC simulation carried out by Nicholas T. Liesen on soft spheres provided a motivation for this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix. Surface tension

The equilibrium density profile of a vapour-liquid interface (at the phase coexistence) is determined by minimising the grand potential $\mathrm{\Omega }[\rho ]$, which is a functional of the density profile. Focusing on the unit area of the interface that runs perpendicular to the z-axis, (A1) $\beta \mathrm{\Omega }[\rho ]={\int }_{-h/2}^{h/2}[f_0(T,\rho )+f_2(T,\rho ){\dot{\rho }}^2-\beta \mu \rho ]\mathrm{d}z,$(A1) where $\dot{\rho }:=\mathrm{d}\rho /\mathrm{d}z$ and h is the height of the system with z = 0 at its centre. We place z = 0 within the interfacial region and suppose that h is taken sufficiently large so that $z=\pm h/2$ is well within the bulk phases.

As discussed in detail in Ref. [19], we can profitably exploit the analogy with mechanics by regarding z as the time and ρ as the position of a particle subject to the potential energy $-f_0+\beta \mu \rho$. The mass of the particle is position dependent and is given by $2f_2(>0)$.

In this mechanical analogy, Equation (A1(A1) $\beta \mathrm{\Omega }[\rho ]={\int }_{-h/2}^{h/2}[f_0(T,\rho )+f_2(T,\rho ){\dot{\rho }}^2-\beta \mu \rho ]\mathrm{d}z,$(A1) ) is the action integral in which (A2) $L(\rho ,\dot{\rho }):=f_0(T,\rho )+f_2(T,\rho ){\dot{\rho }}^2-\beta \mu \rho$(A2) is the Lagrangian. Lagrange's equation of motion (A3) $\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\rho }}\right)-\frac{\mathrm{\partial }L}{\mathrm{\partial }\rho }=0$(A3) leads to (A4) $\beta {\mu }_0(T,\rho )-{\left(\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }\rho }\right)}_T{\dot{\rho }}^2-2f_2(T,\rho )\ddot{\rho }=\beta \mu ,$(A4) where $\ddot{\rho }:={\mathrm{d}}^2\rho /\mathrm{d}{z}^2$ and (A5) $\beta {\mu }_0(T,\rho ):={\left(\frac{\mathrm{\partial }{f}_0}{\mathrm{\partial }\rho }\right)}_T.$(A5) The non-linear second order ordinary differential equation, Equation (A4(A4) $\beta {\mu }_0(T,\rho )-{\left(\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }\rho }\right)}_T{\dot{\rho }}^2-2f_2(T,\rho )\ddot{\rho }=\beta \mu ,$(A4) ), would have to be solved numerically under the boundary conditions that $\underset{z\to \mathrm{\infty }}{lim}\rho ={\rho }_{\mathrm{v}}$ and $\underset{z\to -\mathrm{\infty }}{lim}\rho ={\rho }_{\mathrm{l}}$. This proves to be a rather non-trivial task due to (1) the existence of the additional solutions $\rho (z)\equiv {\rho }_{\mathrm{v}}$ and $\rho (z)\equiv {\rho }_{\mathrm{l}}$ besides the one we seek, (2) an extreme sensitivity of the numerical solution to the initial value of $\dot{\rho }$, say at z = 0 where we may set $\rho =({\rho }_{\mathrm{v}}+{\rho }_{\mathrm{l}})/2$, and (3) the failure of the most of the differential equation solvers to conserve a constant of motion exactly even in a simpler case in which $f_2$ is independent of ρ.

Fortunately, the mechanical analogy mentioned above reveals the existence of the first integral of Equation (A4(A4) $\beta {\mu }_0(T,\rho )-{\left(\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }\rho }\right)}_T{\dot{\rho }}^2-2f_2(T,\rho )\ddot{\rho }=\beta \mu ,$(A4) ). In fact, since L does not depend explicitly on z, the energy function (A6) $e:=\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\rho }}\dot{\rho }-L=f_2{\dot{\rho }}^2-f_0+\beta \mu \rho$(A6) is a constant of motion [31], and hence is independent of z. Evaluating the right-hand-side of Equation (A6(A6) $e:=\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\rho }}\dot{\rho }-L=f_2{\dot{\rho }}^2-f_0+\beta \mu \rho$(A6) ) for either of the bulk phases at coexistence, in which $\dot{\rho }=0$, we find that $e=\beta p$. Recalling the double-tangent construction, by means of which one can find the coexisting densities [31], we see that $f_0+\beta p-\beta \mu \rho \ge 0$. Thus, solving Equation (A6(A6) $e:=\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\rho }}\dot{\rho }-L=f_2{\dot{\rho }}^2-f_0+\beta \mu \rho$(A6) ) for $\dot{\rho }$ and arbitrarily choosing the negative root so that the liquid phase is on the left of the vapour phase, we arrive at [19] (A7) $\dot{\rho }=-\sqrt{\frac{1}{f_2(T,\rho )}[f_0(T,\rho )+\beta p-\beta \mu \rho ]}.$(A7) It is relatively easy to solve this first order ordinary differential equation numerically using, for example, the fourth order Runge-Kutta method under the boundary condition that $\rho =({\rho }_{\mathrm{v}}+{\rho }_{\mathrm{l}})/2$ at z = 0 and integrating in both positive (vapour phase) and negative (liquid phase) directions.

The surface tension is the excess per unit area of the interface of the grand potential over that of the hypothetical system [31–33]. The latter consists of bulk vapour and liquid phases in contact with each other across a sharp dividing surface parallel to the interface, but each behaving as if it is a piece of a bulk homogeneous phase. This leads to the following expression for γ: (A8) $\beta \gamma (T)={\int }_{-h/2}^{h/2}[f_0(T,\rho )+f_2(T,\rho ){\dot{\rho }}^2+\beta p-\beta \mu \rho ]\mathrm{d}z.$(A8) Using Equation (A6(A6) $e:=\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\rho }}\dot{\rho }-L=f_2{\dot{\rho }}^2-f_0+\beta \mu \rho$(A6) ), we can rewrite this equation as (A9) $\beta \gamma (T)=2{\int }_{-h/2}^{h/2}[f_0(T,\rho )+\beta p-\beta \mu \rho ]\mathrm{d}z.$(A9) Another application of Equation (A7(A7) $\dot{\rho }=-\sqrt{\frac{1}{f_2(T,\rho )}[f_0(T,\rho )+\beta p-\beta \mu \rho ]}.$(A7) ) results in Equation (23(23) $\beta \gamma (T)=2{\int }_{{\rho }_{\mathrm{v}}}^{{\rho }_{\mathrm{l}}}\sqrt{f_2(T,\rho )[f_0(T,\rho )+\beta p-\beta \mu \rho ]}\mathrm{d}\rho ,$(23) ). Clearly, it is unnecessary to find $\rho (z)$ if we are interested in γ only.