Andersen-Weeks-Chandler perturbation theory and one-component sticky-hard-spheres

TL;DR

This paper applies second order Andersen-Weeks-Chandler perturbation theory to sticky-hard-spheres fluids and compares its accuracy with other theoretical approaches and Monte Carlo simulations.

Contribution

It demonstrates the effectiveness of second order Andersen-Weeks-Chandler perturbation theory for modeling sticky-hard-spheres fluids.

Findings

Second order perturbation theory aligns well with Monte Carlo results.

Comparison shows advantages over mean spherical and Percus-Yevick approximations.

Provides insights into fluid behavior using advanced perturbation methods.

Abstract

We apply second order Andersen-Weeks-Chandler perturbation theory to the one-component sticky-hard-spheres fluid. We compare the results with the mean spherical approximation, the Percus-Yevick approximation, two generalized Percus-Yevick approximations, and the Monte Carlo simulations.

Equations

$$\beta (A[e]-A[e_0])$$

$$\beta [\Delta A{]}_{(1)}$$

$$\beta [\Delta A{]}_{(2)}$$

$$e(r)=\theta (r-\sigma )+\frac{\sigma }{12\tau }\delta (r-\sigma ),$$

$$e_0(r)=\theta (r-\sigma ),$$

$$\Delta e(r)=\frac{\sigma }{12\tau }\delta (r-\sigma ).$$

$$g_0(1,\dots ,n)\approx \prod _{i<j}^n{g}_0(i,j).$$

$$\beta \frac{[\Delta A{]}_{(1)}}{\bar{N}}$$

$$\beta \frac{[\Delta A{]}_{(2)}}{\bar{N}}$$

$$\frac{I_2}{\rho }=\frac{1}{2\rho }\frac{1}{V}\int d1d2\frac{{\rho }_0(1,2)}{e_0(1,2)}\Delta e(1,2)=\frac{1}{\tau }(\eta {\bar{y}}_0),$$

$$\frac{I_3}{\rho }$$

$$\frac{I_4}{\rho }$$

$$J_3(1,2,3)$$

$$J_4(1,2,3,4)$$

$$\frac{I_3}{\rho }$$

$$\frac{I_4}{\rho }$$

$$g_1$$

$$h_2^a$$

$$h_2^b$$

$$h_3$$

$$h_4$$

$$\cos \delta =\cos {\theta }_1\cos {\theta }_2+\sin {\theta }_1\sin {\theta }_2\cos \varphi .$$

$$a^2$$

$${\bar{y}}_0$$

$$H(z)$$

$$j_0(z)$$

$$\rho {\tilde{c}}_0(z/\sigma )\approx -24\eta {\int }_0^{1}ds{s}^2{j}_0(sz)(\alpha +\beta s+\gamma s^3),$$

$$\alpha$$

$$\beta$$

$$\gamma$$

$$a$$

$${\bar{y}}_0$$

$$\mathcal{C}\bar{N}=\frac{1}{2\bar{N}}\left({\rho }^2\frac{{\chi }_T^0}{{\chi }_T^{id}}\right){\left\{\frac{\partial }{\partial \rho }\beta [\Delta A{]}_{(1)}\right\}}^2,$$

$$\beta \frac{[\Delta A{]}_{(2)}^{\prime }}{\bar{N}}$$

$$\frac{I_3^{\prime }}{\rho }$$

$$\frac{I_4^{\prime }}{\rho }$$

$$\frac{\partial {\rho }_0(1,2)}{\partial \rho }=\frac{1}{\rho }\frac{{\chi }_T^{id}}{{\chi }_T^0}\left\{2{\rho }_0(1,2)+\int d3[{\rho }_0(1,2,3)-\rho {\rho }_0(1,2)]\right\},$$

$$\beta \frac{[\Delta A{]}_{(2)}^{\prime \prime }}{\bar{N}}$$

$$\beta P{\sigma }^3=\beta P_{HS}{\sigma }^3+\frac{6}{\pi }{\eta }^2\frac{\partial \beta \Delta A/N}{\partial \eta },$$

Full text

Andersen-Weeks-Chandler thermodynamic perturbation theory, sticky-hard-spheres, colloidal suspension, mean spherical approximation, Percus-Yevick approximation, generalized Percus-Yevick approximation, Monte Carlo simulation

I Introduction

The sticky-hard-sphere (SHS) model introduced by R. J. Baxter in 1968 Baxter (1968) plays an important role in soft matter offering a description of a sterically stabilized colloidal suspension Fantoni, Gazzillo, and Giacometti (2005); Gazzillo et al. (2006); Gazzillo, Fantoni, and Giacometti (2006); Fantoni et al. (2006, 2007); Gazzillo, Fantoni, and Giacometti (2008, 2009).

In this work we apply Andersen-Weeks-Chandler (AWC) thermodynamic-perturbation-theory (TPT) Hansen and McDonald (1986) to treat the SHS three-dimensional fluid and we compare the results for the equation of state of our calculation with the ones for the mean-spherical-approximation (MSA) Hansen and McDonald (1986), for the Percus-Yevick (PY) approximation Hansen and McDonald (1986), for two generalized-Percus-Yevick (GPY) approximations (C0 and C1 in Ref. Gazzillo and Giacometti (2004)), and for the Monte Carlo simulations of Miller and Frenkel Miller and Frenkel (2003).

We are then able to show how the TPT breaks down at low reduced temperature and high density. Our analysis gives a reference benchmark for the behavior of the SHS system when treated with the AWC TPT scheme.

Our analysis also clarifies the role played by the reducible Mayer diagrams in the second order AWC TPT.

The work is organized as follows. In Section II we introduce the AWC TPT scheme, in Section III we define the SHS fluid model, in Sections IV we outline our calculation of the AWC TPT for the SHS fluid, in Section V we clarify the role played by the reducible integrals, in Section VI we discuss some technical details regarding our Monte Carlo calculation of the various order terms of the TPT, in Section VII we present our results, and Section VIII is for our conclusive discussion.

II The Andersen-Weeks-Chandler thermodynamic perturbation scheme

Following AWC perturbation theory Andersen, Weeks, and Chandler (1971) we consider the Helmholtz free energy $A$ as a functional of the Boltzmann factor $e(1,2)=\exp[-\beta\phi(1,2)]$ ($\phi(1,2)$ being the pair interaction potential of the fluid under exam) and expand it in a Taylor series around the Boltzmann factor, $e_{0}(1,2)$, of a given reference system. Working in the grand-canonical ensemble we obtain the following perturbative expansion in $\Delta e=e-e_{0}$

[TABLE]

where $\beta=1/(k_{B}T)$ (with $k_{B}$ Boltzmann constant and $T$ absolute temperature), $\bar{N}$ average number of particles, $\rho=\bar{N}/V$ (with $V$ volume of the system), $\chi_{T}^{id}=\beta/\rho$ isothermal compressibility of the ideal gas, $\chi_{T}^{0}$ isothermal compressibility of the reference system, $\rho_{0}(1,\ldots,n)$ the grand-canonical ensemble $n-$body correlation function of the reference system, and in the last term of Eq. (3) the density derivative is taken at constant temperature, volume, and chemical potential. In order to derive these expressions one can adapt the details found in Appendix D of Hansen and McDonald book Hansen and McDonald (1986) where their expression (6.2.14) is found. It is then an easy task to pass from their expansion in terms of the pair-potential variation to our expansion in terms of the Boltzmann factor variation.

III One-component sticky-hard-spheres

For the Baxter Baxter (1968) one-component sticky-hard-spheres (SHS) model one has

[TABLE]

where $\sigma$ is the spheres diameter, $\tau$ the reduced temperature, $\theta$ is the Heaviside step function, and $\delta$ the Dirac delta function.

Choosing as reference system the hard-spheres (HS) model one has

[TABLE]

so that

[TABLE]

So one sees that AWC expansion (1) reduces to an expansion in powers of $1/\tau$.

IV Calculation

Before expression (3) can be used some approximation must be introduced for the three- and four-body distribution functions. The most widely used approximation is Kirkwood superposition approximation Kirkwood (1935). This has previously successfully applied to the second order thermodynamic perturbation study of the square well potential by Henderson and Barker Henderson and Barker (1971).

Using the Kirkwood superposition approximation (KSA) Kirkwood (1935) one can express the $n-$body correlation functions $\rho_{0}(1,\ldots,n)=\rho^{n}g_{0}(1,\ldots,n)$ in terms of pair distribution functions according to

[TABLE]

The idea is to use for the pair distribution function of the reference HS system the analytic solution of the Ornstein-Zernike equation with the Percus-Yevick closure.

The first two terms in the perturbative expansion (1) reduce to

[TABLE]

where

[TABLE]

where $\eta=\frac{\pi}{6}\rho\sigma^{3}$ is the hard sphere packing fraction, $y_{0}(1,2)=g_{0}(1,2)/e_{0}(1,2)$ is the cavity function of the reference system and $\bar{y}_{0}=y_{0}(|{\mathbf{r}}_{1}-{\mathbf{r}}_{2}|=\sigma)$. Upon using KSA one finds,

[TABLE]

where we have introduced

[TABLE]

where $h_{0}(1,2)=g_{0}(1,2)-1$ is the total correlation function of the reference system. Note that the first term in $J_{3}$ and the first and second terms in $J_{4}$ give rise to reducible integrals (i.e. integrals that can be reduced into products of simpler integrals).

It is convenient to perform the calculation of $I_{3}$ and $I_{4}$ in reciprocal space, to get,

[TABLE]

and

[TABLE]

where in the integrand of $h_{4}$

[TABLE]

In all these expressions we have introduced the following notation

[TABLE]

where $g_{0}(r),y_{0}(r),\tilde{h}_{0}(k),\tilde{c}_{0}(k)$ are respectively the hard spheres radial distribution function, cavity function, the Fourier transform of the total correlation function and the Fourier transform of the direct correlation function, and $j_{0}$ is the zeroth order spherical Bessel function of the first kind.

Finally the Fourier transform of the HS direct correlation function calculated through the Percus-Yevick closure is given by Ashcroft and Lekner (1966)

[TABLE]

where

[TABLE]

and it is easily verified that under such approximation one has

[TABLE]

V Neglecting reducible integrals

It has been observed by Henderson and Barker Henderson and Barker (1971) that the role of the last term in Eq. (3)

[TABLE]

is to cancel in the second order term of the perturbative expansion, $[\Delta A]_{(2)}$, all reducible integrals appearing in $I_{3}$ and $I_{4}$. So that the final expression for the second order term of expansion (1) would be (exactly the expression found in Andersen, Weeks, and Chandler (1971))

[TABLE]

where

[TABLE]

Alternatively one may use the sum rule

[TABLE]

to rewrite ${\cal C}$ [Eq. (33)] in terms of two and three body correlation functions and upon using the superposition approximation one finds

[TABLE]

VI Technical details

The five integrals (17)-(21) where all calculated using Monte Carlo technique Kalos and Whitlock (1986) averaging the various integrands on $10^{6}$ randomly sampled points. Since all of those integrals are improper (extending up to infinity in the $z$ variables) it was necessary to split each integration on the $z$ variables into an integral over $[0,1]$ plus an integral over $[1,\infty]$. This latter integral was then reduced through a change of variable $z\to 1/z$ into an integral over $[0,1]$.

The errors on the estimate of a given integral was calculated so that the true value of the integral would lie $99.7\%$ of the time within the estimate plus or minus the error.

VII Results

Figs. 1-4 show the results for $\beta\Delta A/N$ as a function of $\eta$. Amongst the three expressions used: (9), (34), and (38)), the more accurate is $[\Delta A]_{(2)}^{\prime}$, the one suggested in Andersen, Weeks, and Chandler (1971) and it falls on the PY approximation for big $\tau$ and small $\eta$. At high $\eta$ the error bars become more relevant.

Figs. 8-5 show the results for

[TABLE]

as a function of $\eta$, where for the pressure of the HS reference system we chose the PY result from the compressibility route, i.e.

[TABLE]

The second order AWC TPT is taken from the (34 calculation.

VIII Discussion

Our first calculation, the one using $[\Delta A]_{(2)}$ (see Eq. (9)) is certainly not correct because we are using the KSA only on the integrands of the first two integrals of Eq. (3) calculating the last term exactly; this certainly leads to an inconsistency in the use of KSA.

Our third calculation, the one using $[\Delta A]_{(2)}^{\prime\prime}$ (see Eq. (38)) is also not correct. This can be understood as follows. It is well known that KSA fails to satisfy the sum rule (37). Using KSA in the left hand side of Eq. (37) one finds

[TABLE]

where

[TABLE]

and we used the compressibility sum rule,

[TABLE]

Eq. (41) can be also rewritten as,

[TABLE]

This approximation is certainly valid in the limit of small densities when $\chi_{T}^{0}\to\chi_{T}^{id}$ and $h_{0}\to e_{0}-1=f_{0}$ ($f_{0}$ being the Mayer function of the reference system), after all the KSA becomes exact in such limit (as the potential of mean force tends to the pair interaction potential). Otherwise the correction term $\gamma_{3}/(\rho^{2}g_{0})$ would be of order $\rho$ as $\rho\to 0$ (see the Appendix). So that the exact expression for the density derivative of the two body correlation function would be

[TABLE]

where $\gamma_{i}=O(\rho^{i})$ as $\rho\to 0$. It is then clear that in calculating the square

[TABLE]

in the ${\cal C}$ term, the term stemming from

[TABLE]

which gives rise to the last term in Eq. (38), will be of the same leading order ($\rho^{4}$) as the one coming from

[TABLE]

in the small density limit. But since in KSA this last term is neglected, in order to be consistent (up to orders $\rho^{3}$ in the small density limit) one needs to neglect also the term of Eq. (48). Moreover it can be easily verified that the two terms coming from $\gamma_{1}$ times $\gamma_{1}$ cancel the first reducible integral in $I_{3}$ and the first reducible integral in $I_{4}$ whereas the term coming from $\gamma_{1}$ times $\gamma_{2}$ cancels the second reducible integral in $I_{4}$. So that Eq. (34) (the original AWC expression) for the second order perturbative term in the AWC theory, is recovered.

The correct second order AWC calculation, $[\Delta A]_{(2)}^{\prime}$ (see Eq. (34)) shows that the TPT breaks down at small reduced temperatures $\tau$ and large packing fractions $\eta$, as expected.

Appendix A Correction to approximation (45)

One can understand that Eq. (45) is not an exact relation by comparing the small density expansion of the left and right hand side. For the left hand side we have Hansen and McDonald (1986)

[TABLE]

where in the Mayer graphs the filled circles are field points of weight $1$ and connecting bonds are Mayer functions of the reference system $f_{0}$. And using

[TABLE]

in the right hand side one finds,

[TABLE]

So that the correction term is of order $\rho$, namely,

[TABLE]

The correct small density expansion for the density derivative of the two body correlation function is

[TABLE]

where the first term neglected in KSA is $\rho^{2}\alpha_{1}^{\prime}=O(\rho^{3})$.