---
abstract: 'It has recently been shown that a free energy for Baxter’s sticky hard sphere fluid is uniquely defined within the framework of fundamental measure theory (FMT) for the inhomogeneous hard sphere fluid, provided that it obeys scaled-particle theory and the Percus-Yevick (PY) result for the direct correlation function \[Hansen-Goos and Wettlaufer, J. Chem.  Phys. [** 134**]{}, 014506 (2011)\].  Here, combining weighted densities from common versions of FMT with a new vectorial weighted density, we derive a regularization of the divergences of the associated strongly confined limit.  Moreover, the simple free energy that emerges is exact in the zero-dimensional limit, leaves the underlying equation of state unaffected, and yields a direct correlation function distinct from the PY expression.  Comparison with simulation data for both the bulk pair correlation function and the density profiles in confinement shows that the new theory is significantly more accurate than the PY-based results.  Finally, the resulting free energy is applicable to a glass of adhesive hard spheres.' 
author:
- 'Hendrik Hansen-Goos'
- 'Mark A. Miller'
- 'J.S. Wettlaufer'
title: 'Density functional theory for Baxter’s sticky hard spheres in confinement'
---

Colloidal suspensions are known to exhibit a variety of collective phenomena such as crystallization, gel formation, glass transition and percolation [@BuRuPi07; @Luetal08; @EbWaCa11].  Unlike in molecular systems these properties can be measured not only in terms of bulk sampling such as x-ray diffraction, heat capacity measurements, or mechanical properties, but also on a microscopic level thanks to optical single-particle observation techniques such as confocal microscopy.  An ideal theoretical tool to study the structural properties down to the particle size is classical density functional theory (DFT) [@Ev79], an excellent example application being the purely entropic hard-sphere system.  Indeed, within the framework of fundamental measure theory (FMT) [@Ro89; @KiRo90; @KiRo93], for arbitrary external potentials both the inhomogeneous liquid and the crystalline phases are described by a simple and accurate free energy.  However, a more realistic model of colloidal suspensions must account for attractive interactions.  Such a model is provided by Baxter’s sticky hard spheres (SHS) which interact via an infinitely narrow and deep attractive well, the limit being taken such that the second virial coefficient remains finite [@Ba68].  The model is particularly appealing because it is soluble analytically within the Percus Yevick (PY) approximation.  The PY results obtained by Baxter provide a framework to interpret experimental phase diagrams of colloidal suspensions, including the liquid-liquid transition, crystallization, percolation and, to some extent, the glass transition, revealing the existence of attractive and repulsive glasses [@Phetal02].  Quantitative interpretation of single-particle observations in colloidal suspensions rely on a theoretical approach capable of resolving the microscopic structure of Baxter’s SHS as accurately as does FMT for hard spheres.  Previous DFT approaches relied on Taylor expansions of the free energy around the homogeneous bulk fluid which use the PY result for the direct correlation function while crudely approximating higher order terms.  Recently, using a set of weighted densities from FMT for the hard-sphere fluid [@HaWe11], a first non-perturbative DFT was derived.  When applied to density profiles near a planar wall, or the SHS confined between two parallel planar hard walls, the theory was shown to be a significant improvement over previous approaches.  However, based on experience with Rosenfeld’s original FMT for the hard-sphere fluid, it is expected that the FMT for the SHS would exhibit an unphysical divergence for strongly confined fluids, the most extreme test being the zero-dimensional (0D) limit. 

In this Letter, we recast the original FMT for the SHS, which is based on Kierlik’s and Rosinberg’s scalar weighted densities [@KiRo90; @KiRo93], by combining Rosenfeld’s vectorial weighted density [@Ro89] with a new vectorial weighted density.  In consequence, we identify the origin of the divergence and regularize it in the spirit of Schmidt [*et al.*]{} ’s approach to Rosenfeld’s FMT for the hard-sphere system [@Roetal97].  Rather than deriving a cumbersome expression by conserving the underlying PY correlation function, we seek a [*simple*]{} density functional.  When the stickiness of the particles is strong, the modified theory (exact in the 0D limit) provides a better description of the pair-correlation function than does the PY-based FMT, while the accurate results at moderate stickiness are left virtually unchanged.  For the SHS confined in a spherical cavity, the new theory is significantly better than its predecessor.  Our free energy can be used to study a wide range of the structural properties of inhomogeneous colloidal suspensions with short-ranged attractions.  In particular, this is the first free energy applicable to the glassy state where confinement due to caging is ubiquitous.  In fact, it may overcome previous theoretical pathologies associated with Baxter’s adhesive hard-sphere treatment of the glass transition.  These problems were based on using the PY direct correlation function alone [@Foetal00], rather than the complete free energy provided here. 

Consider a fluid of monodisperse hard spheres with radius $R$. We define the locally averaged packing fraction $n_3({{\bf\stringr}}) = \int{\mathrm{d}}{{\bf\stringr}}'\rho({{\bf\stringr}}')\Theta(R-|{{\bf\stringr}}-{{\bf\stringr}}'|)$ with weighted densities $n_2({{\bf\stringr}}) = \frac{\partial}{\partial R} n_3({{\bf\stringr}})$ and ${{\bf\stringn}}_2({{\bf\stringr}}) = -\nabla n_3({{\bf\stringr}})$. Rosenfeld’s excess free energy density (measured in units of $k_B T$) is then $\Phi_{\text{HS}} = \Phi_1 + \Phi_2 + \Phi_3$ with $\Phi_1 =  -  n_2 \ln(1-n_3)/(4\pi R^2)$, $\Phi_2 = (n_2^2 - {{\bf\stringn}}_2\cdot{{\bf\stringn}}_2)/[4\pi R (1-n_3)]$, and $\Phi_3 = (n_2^3-3 n_2 {{\bf\stringn}}_2\cdot{{\bf\stringn}}_2)/[24\pi(1-n_3)^2]$ [@Ro89].  While $\Phi_{\text{HS}}$ has been extremely successful in describing an [*unconfined*]{} fluid, with an inhomogeneity caused, say, by a planar wall, it is problematic when the fluid is highly confined (see [@Ta00] and refs.  therein).  For the most extreme confinement, the density distribution becomes simply a delta-peak, i.e., $\rho^{\text{0D}}(r) = \eta \delta(r)$, where $r = |{{\bf\stringr}}|$, and $\eta$ is the packing fraction ($0\le \eta \le 1$).  The weighted densities for this distribution are $n_3^{\text{0D}}(r) = \eta\Theta(R-r)$, $n_2^{\text{0D}}(r) = \eta\delta(R-r)$, and ${{\bf\stringn}}_2^{\text{0D}}({{\bf\stringr}}) = \eta\delta(R-r) \hat{{{\bf\stringr}}}$, where $\hat{{{\bf\stringr}}}={{\bf\stringr}}/r$. The integration of $\Phi_1$ can be readily performed yielding $\int{\mathrm{d}}{{\bf\stringr}} \Phi_1 = \eta + (1-\eta) \ln(1-\eta)$ which is the exact result for a 0D cavity holding an average number of particles $\eta$ [@Roetal97].  Examining $\Phi_{\text{HS}}$ further for such a cavity shows that $\Phi_2$ vanishes (as it should) while $\Phi_3$ exhibits a strong negative divergence.  This is due to the factor $3$ in front of the vector term in $\Phi_3$ which would have to be unity for its scalar and vectorial terms to cancel.  However, the resulting expression for the excess free energy does not yield the PY direct correlation function, thereby destroying the excellent experimental agreement of Rosenfeld’s $\Phi_{\text{HS}}$ for the HS fluid.  Several routes have been devised to regularize $\Phi_3$. A particularly appealing approach consists of introducing a new weighted density $\bar{n}_2 = n_2-({{\bf\stringn}}_2\cdot{{\bf\stringn}}_2)/n_2$ that reduces to $n_2$ in the bulk fluid (where the vectorial weighted densities vanish as $n_3$ is spatially uniform), and vanishes in the 0D limit (where ${{\bf\stringn}}_2\cdot{{\bf\stringn}}_2= n_2^2$).  Using $\bar{n}_2$, a regularized excess free energy density proposed by Schmidt [*et al.*]{} [ @Roetal97] is $$\bar{\Phi}_{\text{HS}} = - \frac{n_2 \ln(1-n_3)}{4\pi R^2} + \frac{n_2 \bar{n}_2}{4\pi R(1-n_3)} + \frac{\bar{n}_2^3}{24\pi(1-n_3)^2} \, .$$

From the properties of $\bar{n}_2$ it is immediately clear that $\bar{\Phi}_{\text{HS}}$ is exact for a 0D cavity while the free energy of the bulk fluid is unaffected.  However, the third term differs by order $[({{\bf\stringn}}_2\cdot {{\bf\stringn}}_2)/n_2^2]^2$ from that of $\Phi_{\text{HS}}$, which is very small for moderately inhomogeneous fluids.  In fact, it affects only the direct correlation functions of order 4 and higher whereas the PY direct correlation function (which is second order) is preserved.  Both unconfined and confined fluids are extremely well described by $\bar{\Phi}_{\text{HS}}$; accurately predicting solid-liquid coexistence in hard sphere systems whereas $\Phi_{\text{HS}}$ dramatically overstabilizes the crystal due to the negative divergence for strongly peaked density distributions.  However, Tarazona’s modification of $\Phi_3$, which uses a tensorial weighted density [@Ta00], performs better with regard to the width of the density peaks in the crystal. 

How does the picture change as we introduce Baxter’s attractive surface interaction into the system?  Rosenfeld’s set of weighted densities (i.e. $n_3$, $n_2$, and ${{\bf\stringn}}_2$) is known to be insufficient for constructing any reasonable free energy density for the SHS fluid because they cannot yield the $\delta$-function within the direct correlation function [@HaWe11].  However, when $n_3$ and $n_2$ are supplemented by the weighted density $n_1({{\bf\stringr}}) = \frac{1}{8\pi}\frac{\partial}{\partial R} n_2({{\bf\stringr}})$, the set is sufficient to construct a free energy density $\Phi_{\text{SHS}} = \Phi_{\text{HS}} + \Phi_{\text{S}}$ for Baxter’s SHS fluid where $\Phi_{\text{S}}=n_1\phi_1/R + n_2\phi_2/(2\pi R^2)$ is uniquely defined by (a) requiring the density functional to yield the PY direct correlation function and (b) imposing consistency with scaled particle theory [@HaWe11].  The dimensionless coefficients $\phi_1$ and $\phi_2$ are functions of $x = \frac{R n_2}{1-n_3}$, which in the bulk fluid is related to the packing fraction $\eta$ via $x = \frac{3\eta}{1-\eta}$. They are obtained by solving $\phi_1'(x) = - 2 \tilde{y}_\sigma(x)/x$ and $[x^2 \phi_2'(x)]' = \tilde{y}_\sigma(x)^2/2 - x \tilde{y}_\sigma(x) + x \tilde{y}'_\sigma(x)$, where the integration constants are chosen such that $\phi_1$ and $\phi_2$ vanish for $x\to 0$. Here $\tilde{y}_{\sigma}=\eta y_{\sigma}/\tau$, where $y_{\sigma}$ is the cavity function at contact, and is obtained as the smaller of the two solutions to $\tilde{y}_{\sigma}^2-(4x+12\tau)\tilde{y}_{\sigma}+2x(2+x)=0$.

In order to regularize the divergences in the 0D limit we rewrite $\Phi_{\text{S}}$ using vectorial weighted densities so that the resulting functional $\bar{\Phi}_{\text{S}}$ vanishes.  This guarantees that $\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}}$ is exact for a cavity that can hold only one particle at a time, and thus the corresponding free energy is unaffected by the sticky interaction.  We use the identity $n_1({{\bf\stringr}}) = n_2({{\bf\stringr}})/(4\pi R) + (\nabla\cdot\nabla)n_3({{\bf\stringr}})/(8\pi)$ [@KiRo93] in $\Phi_{\text{S}}$, and for reasons that will become clear below, we split the term $n_1\phi_1/R$ in half and substitute for $n_1$ in only one of the terms.  Moreover, because $\Phi_{\text{S}}$ is always integrated on $\mathbb{R}^3$, we can apply Green’s first identity and neglect the boundary integral.  The result is $\Phi_{\text{S}} = \Phi_{\text{S1}}+\Phi_{\text{S2}}$, where $\Phi_{\text{S1}}= \left(\phi_{12} n_1 n_2-  \phi_{12}^{\text{v}}{{\bf\stringn}}_1\cdot {{\bf\stringn}}_2\right)/(1-n_3)$ and $\Phi_{\text{S2}}= \left[\phi_{222} n_2^3-2 \phi_{222}^{\text{v}} n_2 ({{\bf\stringn}}_2\cdot {{\bf\stringn}}_2)\right]/[4\pi(1-n_3)^2]$. The coefficients $\phi_{12} = \phi_1/(2x)$, $\phi_{12}^{\text{v}} = \phi_1'/2$, $\phi_{222} = (\phi_1/2 + 2 \phi_2)/x^2$, and $\phi_{222}^{\text{v}} = \phi_1'/8$ are functions of $x = \frac{R n_2}{1-n_3}$, and the new vectorial weighted density is ${{\bf\stringn}}_1({{\bf\stringr}}) = -\nabla n_2({{\bf\stringr}})/(8\pi)$.

Replacing $n_1$ only in half of the split term $n_1\phi_1/R$ has two important consequences.  Firstly, it insures that the function $\phi_{222}$ remains finite in the limit $x\to 0$. Secondly, it leads to a symmetry between $\phi_{12}$ and $\phi_{12}^{\text{v}}$ that is reflected in their equality in the limits $x \to 0$ and $x\to\infty$, where $\phi_{12}(0)=\phi_{12}^{\text{v}}(0)=-\frac{1}{3\tau}$ and $\phi_{12}(\infty)=\phi_{12}^{\text{v}}(\infty)=\sqrt{2}-2$. However, comparing $\phi_{12}$ with $\phi_{12}^{\text{v}}$ for any $x$ reveals their differences (Fig.  \[fig\_phis\]).  When we replace $\phi_{12}^{\text{v}}$ with $\phi_{12}$, $\Phi_{\text{S1}}$ is unaffected in the bulk fluid.  Above we calculated the weighted densities $n_3^{\text{0D}}$, $n_2^{\text{0D}}$ and ${{\bf\stringn}}_2^{\text{0D}}$ in the 0D limit.  These are now supplemented by $n_1^{\text{0D}}(r) = \eta\delta'(R-r)/(8\pi)$, and ${{\bf\stringn}}_1^{\text{0D}}({{\bf\stringr}}) = \eta\delta'(R-r) \hat{{{\bf\stringr}}}/(8\pi)$, from which we have $n_1^{\text{0D}}n_2^{\text{0D}}-{{\bf\stringn}}_1^{\text{0D}}\cdot {{\bf\stringn}}_2^{\text{0D}}\equiv 0$. Hence, this modification regularizes $\Phi_{\text{S1}}$. Note, that because $\phi_{12}^{\text{v}}-\phi_{12}$ vanishes as $x^{-1}\ln x$ for $x\to\infty$ then even if $\phi_{12}^{\text{v}}$ and $\phi_{12}$ take the same values in this limit it is insufficient to insure that $\Phi_{\text{S1}}$ vanishes in the 0D limit, where large values of $x$ occur due to the delta-function in $n_2^{\text{0D}}$. Furthermore, larger powers of the vectorial weighted densities could be used to modify the terms $\phi_{12} n_1 n_2/(1-n_3)$ and $\phi_{12}^{\text{v}}{{\bf\stringn}}_1\cdot {{\bf\stringn}}_2/(1-n_3)$ individually while maintaining for example the PY direct correlation function which underlies $\Phi_{\text{S}}$. However, the symmetry of $\phi_{12}$ and $\phi_{12}^{\text{v}}$ motivates the simple remedy introduced here.  The value of giving up Baxter’s PY result for the direct correlation function is assessed in terms of the performance of the modified functional in comparison to simulations. 

The regularization of $\Phi_{\text{S2}}$ requires two steps.  We replace $\phi_{222}^{\text{v}}$ by $\phi_{222}$ and note that even for asymptotic values of the argument $x$, $\phi_{222}(0) = -\frac{1}{12\tau}+\frac{1}{648\tau^3}$ while $\phi_{222}^{\text{v}}(0)=-\frac{1}{12\tau}$; and $\phi_{222}(\infty) = -\frac{1}{3}(\sqrt{2}-1)\simeq -0.138$ while $\phi_{222}^{\text{v}}(\infty)=-\frac{1}{4}(2-\sqrt{2})=-0.146$, these function differ (see Fig.  \[fig\_phis\]).  Therefore, our substitution is justified by the fact that it becomes exact in the limit of low density ($x\to0$) [*and*]{} weak “stickiness” ($\tau\to\infty$).  Like $\Phi_{\text{S1}}$, $\Phi_{\text{S2}}$ in the bulk fluid remains unaffected by the substitution. 

After replacing $\phi_{222}^{\text{v}}$ by $\phi_{222}$, $\Phi_{\text{S2}}$ is still divergent.  This is a consequence of the factor $2$ in front of the vectorial term, which comes from matching the $x\to 0$ and $\tau\to\infty$ limits of $\phi_{222}$ and $\phi_{222}^{\text{v}}$. In analogy to $\bar{\Phi}_{\text{HS}}$ where the term $n_2^3-3n_2({{\bf\stringn}}_2\cdot{{\bf\stringn}}_2)$ is regularized by replacing it with $\bar{n}_2^3$ we now regularize $n_2^3-2n_2({{\bf\stringn}}_2\cdot{{\bf\stringn}}_2)$ by substituting $n_2\bar{n}_2^2$. Introducing the modified weighted density $\bar{n}_1 = n_1- ({{\bf\stringn}}_1\cdot{{\bf\stringn}}_2)/n_2$ (which equals $n_1$ in the bulk fluid and vanishes in the 0D limit) we can write down the regularized version of $\Phi_{\text{S}}$ as $$\label{eq_PhiSreg}
 \bar{\Phi}_{\text{S}} =\frac{\phi_{12}\bar{n}_1 n_2}{1-n_3}
   + \frac{\phi_{222}n_2 \bar{n}_2^2}{4\pi(1-n_3)^2} \, ,$$ where $\phi_{12}$ and $\phi_{222}$ are functions of $x = \frac{R n_2}{1-n_3}$.

As a first test of this regularized free energy density $\bar{\Phi}_{\text{S}}$ we calculate the pair correlation function $g(r)$ by minimizing the density functional $\Omega = \mathcal{F}_{\text{id}} + \int{\mathrm{d}}{{\bf\stringr}} [\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}}+\rho({{\bf\stringr}})(V_{\text{ext}}({{\bf\stringr}})-\mu)]$, where $\mathcal{F}_{\text{id}}$ is the functional for the ideal gas, $\mu$ is the chemical potential and $V_{\text{ext}}$ is the external potential, chosen so that it represents a particle of the SHS fluid at position $r=0$. The resulting density profile $\rho(r)$ gives $g(r) = \rho(r)/\rho_{\text{b}}$, where $\rho_{\text{b}}$ is the bulk fluid density.  In Fig.  \[fig\_gr\] we show the case of strong adhesion, $\tau = 0.13$, and compare the FMT results with MC simulations obtained using the algorithm described in [@MiFr04].  The regularized $\bar{\Phi}_{\text{S}}$ is superior to the PY-based $\Phi_{\text{S}}$, which is remarkable because $\bar{\Phi}_{\text{S}}$ has been optimized for the 0D limit and not at all for $g(r)$. Other, larger $\tau$, MC simulations for $g(r)$ (see [@MiFr04b]) are equally well described by the two FMTs, which give very similar results for $\tau>0.2$.

Alternatively, an analytical expression for the direct correlation function $c(r)$ can be obtained from the excess free energy $\mathcal{F}_{\text{ex}} = \int {\mathrm{d}}{{\bf\stringr}}(\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}})$ through $$\label{eq_DCF}
  c(r) = - \frac{\delta^2\mathcal{F}_{\text{ex}}}{\delta \rho({{\bf\stringr}}) \delta \rho({{\bf\stringr}}')} = - \sum_{\alpha, \beta} \frac{\partial^2(\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}})}{\partial n_{\alpha} \partial n_{\beta} } \omega_{\alpha \beta}(r) \, ,$$ where $\alpha$ and $\beta$ run through $\{ 1, 2, 3, {\bf 1}, {\bf 2}\}$ and $r = |{{\bf\stringr}}-{{\bf\stringr}}'|$. The partial derivatives of $\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}}$ are evaluated in the bulk fluid and derivatives with respect to vectorial quantities ${{\bf\stringn}}_1$ and ${{\bf\stringn}}_2$ are executed formally while the fact that ${{\bf\stringn}}_1=0$ and ${{\bf\stringn}}_2=0$ in bulk ensures that the resulting expressions are scalars.  The symmetric coefficients $\omega_{\alpha \beta}$ are convolutions of the weight functions that are associated with the respective weighted densities.  Therefore $\omega_{\alpha \beta}\equiv 0$ for $r>2R$. For example, $\omega_{33} = \frac{\pi}{12}(4R+r)(2R-r)^2$ is the overlap volume of two spheres with radius $R$ and center-to-center distance $r$. Moreover, $\omega_{23} = \pi R(2R-r)$, $\omega_{13} =-r/8+R/2-R^2/(4r)$, $\omega_{22} =2\pi R^2/r$, $\omega_{12} =R/(4r)+R\delta(2R-r)/8$, $\omega_{\bf 22} = 2\pi R^2/r-\pi r$, and $\omega_{\bf 12} =R/(4r)- R \delta(2R-r)/8$ are required in Eq.  . 

In particular, the coefficient of the $\delta$-function in $c(r)$ is $a=-\frac{R}{4}\left(\frac{2\phi_{12}}{1-\eta}+\frac{3\eta\phi_{12}'}{(1-\eta)^2} \right)$. This result can be used to calculate the percolation line which, under the approximations described in [@ChGl83], is $\frac{R}{6a} = \eta$. As long as $\tau$ is not too large, both the classic PY result and the modified FMT describe the percolation threshold as determined in simulations, the latter being somewhat more reliable in the range of intermediate $\tau$ (Fig.  \[fig\_perc\]). 

Finally, we consider a confinement scenario for which our theory is expected to be the most relevant.  We compare the previous PY-based FMT and the present regularized FMT to canonical MC simulations for a SHS fluid with $\tau=0.13$ and $\tau=0.2$ in a spherical cavity of radius $R_{\text{cav}}=3.7 R$. This system was considered previously [@Goetal97] for a hard-sphere fluid ($\tau\to\infty$).  To compare the FMTs with the simulations we have to transform the former from the grand canonical ensemble to the canonical ensemble which can be achieved up to corrections of order $1/N^3$ by using the scheme of [@Goetal97] that is based on results from [@Saetal96].  Fig.  \[fig\_profiles\] shows that the regularized FMT is clearly superior to the previous PY-based FMT, especially in the center of the spherical cavity.  Interestingly, while the previous FMT underestimates the density in the center when $\tau$ is small, it overestimates the density in the center when applied to hard spheres; large $\tau$ (see [@Goetal97]).  The crossover occurs around $\tau=0.5$ where the regularized and the PY-based FMT yield very similar results.  In the limit of large cavity radius $R_{\text{cav}}$, planar wall density profiles are recovered; comparison (not shown) with existing MC simulations [@JaBr94] shows that the PY-based and the regularized FMT describe the density profiles equally well.  In particular, the theories yield identical contact values $\rho_c$ for the density at a hard wall.  This results from $\bar{\Phi}_{\text{HS}}+\bar{\Phi}_{\text{S}}$ being equal to $\Phi_{\text{HS}}+\Phi_{\text{S}}$ in bulk and hence the underlying pressure being the PY compressibility result $p_{\text{PY}}$ for both theories.  Given that both FMTs obey the contact theorem $\rho_c=\beta p_{\text{PY}}$ the contact values of the density must be identical.  This illustrates the importance of preserving the bulk fluid properties while removing divergences for highly peaked density distributions in order to obtain a density functional theory that is robust and accurate for a broad spectrum of settings. 

[1]{}

S. Buzzaccaro, R. Rusconi, and R. Piazza, Phys.  Rev. Lett. [** 99**]{}, 098301 (2007). 

P.J. Lu, E. Zaccarelli, F. Ciulla, A.B. Schofield, F. Sciortino, and D.A. Weitz, Nature [**453**]{}, 499 (2008). 

A.P.R. Eberle, N.J. Wagner, and R. Casta[ñ]{}eda-Priego, Phys.  Rev. Lett. [** 106**]{}, 105704 (2011). 

R. Evans, Adv.  Phys. [** 28**]{}, 143 (1979). 

Y. Rosenfeld, Phys.  Rev. Lett. [** 63**]{}, 980 (1989). 

E. Kierlik and M.L. Rosinberg, Phys.  Rev. A [**42**]{}, 3382 (1990). 

S. Phan, E. Kierlik, M.L. Rosinberg, B. Bildstein, and G. Kahl, Phys.  Rev. E [**48**]{} 618 (1993). 

R.J. Baxter, J. Chem.  Phys. [** 49**]{}, 2770 (1968). 

K.N. Pham, A.M. Puertas, J. Bergenholtz, S.U. Egelhaaf, A. Moussa[ï]{}d, P.N. Pusey, A.B. Schofield, M.E. Cates, M. Fuchs, and W.C.K. Poon, Science [**296**]{}, 104 (2002). 

H. Hansen-Goos and J.S. Wettlaufer, J. Chem.  Phys. [** 134**]{}, 014506 (2011). 

Y. Rosenfeld, M. Schmidt, H. Löwen, and P. Tarazona, Phys.  Rev. E [**55**]{}, 4245 (1997). 

G. Foffi, E. Zaccarelli, F. Sciortino, P. Tartaglia, and K.A. Dawson, J. Stat.  Phys. [** 100**]{}, 363 (2000). 

P. Tarazona, Phys.  Rev. Lett. [** 84**]{}, 694 (2000). 

M.A. Miller and D. Frenkel, J. Chem.  Phys. [** 121**]{}, 535 (2004). 

M.A. Miller and D. Frenkel, J. Phys.:  Condens.  Matter [**16**]{}, S4901 (2004). 

S.B. Lee, J. Chem.  Phys. [** 114**]{}, 2304 (2001). 

Y.C. Chiew and E.D. Glandt, J. Phys.  A: Math.  Gen. [**16**]{}, 2599 (1983). 

A. Gonz[á]{}lez, J.A. White, F.L. Rom[á]{}n, S. Velasco, and R. Evans, Phy.  Rev. Lett. [** 79**]{}, 2466 (1997). 

J.J. Salacuse, A.R. Denton, and P.A. Egelstaff, Phys.  Rev. E [**53**]{}, 2382 (1996). 

A. Jamnik and D. Bratko, Phys.  Rev. E [**50**]{}, 1151 (1994).