text:
# Phase equilibria and glass transition in colloidal systems with short-ranged attractive interactions. Application to protein crystallization.

## I Introduction

Recently novel results in the statistical mechanics of the fluid state have emerged from a series of studies of the phase diagram of particles with a hard core and short-ranged attractive potential. At first sight this type of interaction potential would have been expected to yield the well known phase diagrams of simple solids, liquids and gases. Here we shall see that a new scenario emerges.
In many practical situations attractive forces that are short-ranged compared to the size of the particles arise because the size of the particles itself is large, whilst the physical forces retain their typical microscopic range. Thus, in order to model large molecules, such as proteins, colloids and nano-particles, one often works in a regime where the ratio of the range of attraction to the size of the repulsive core is small. This crucial issue is emerging in the literature in such areas as protein crystallization , dense colloids , nano-particle assemblies, pre-ceramic particle gelation, latex formation, “buckyballs”  and many others. In some arenas its relevance begins to be recognized, while in others it still remains to be understood.
Atoms or molecules interacting via a potential with a repulsive core and an attractive tail may form gases, liquids and solids as a function of the temperature and density. The attraction is responsible for the liquid-gas transition, while the solid is dominated by the repulsion. In the low-to-medium range of densities, particle entropy and energy enter the free energy and compete to produce liquid (energy-favoured) and gas (entropy-favoured) phases. Thus, beneath the critical point, the total system free energy is optimized if the system splits into two sub-systems, and thereby a sum of two independent free-energies. One of these free energy terms is dominated by the energy (liquid) and the other (gas) by the entropy. A homogeneous, unseparated system, maintained at the average density of the two sub-systems, would not be able to fully optimize either energy or entropy. This, in essence, is the cause of typical first-order phase transitions between a liquid and a gas. In the liquid, as the density increases the remaining diffusive motions of the particles, are reduced and their attendant (configurational) entropy diminishes. Now, the random structure of a liquid is favored provided overall diffusive particle motion is possible. When, at higher densities, diffusion is greatly limited, the free volume is better utilized by making a regular crystalline array, and the entropy of the system is increased by the vibrations of particles within the regular array of a crystal. This is the reason why hard-sphere fluids begin to form a regular crystal at about $`49\%`$ volume fraction, and above the freezing transition at $`55\%`$ have a higher entropy than a liquid-like structure at the same volume fraction. Indeed, for hard spheres, at volume fraction between $`49\%`$ and $`55\%`$ the system again phase separates to achieve the optimal value of entropy: a lower density fluid of $`49\%`$ volume fraction and a well-packed crystal at $`55\%`$. Nevertheless, it is well known that in some experimental studies the system may be unable to access the crystalline state rapidly enough, and the disordered super-cooled liquid structure freezes into a glass. For hard sphere fluids, this occurs at a volume fraction higher than about $`58\%`$. The formation of the glass is a signal that equilibrium quantities such as the free energy may not be sufficient, to describe the behavior of the system, and that long-lived dynamically arrested states may be important under such conditions.

It is possible to relate these qualitative comments to recent understanding acquired from experiments, simulations and theory. Thus, there is a separation of time scales in dense systems, including super-cooled liquids, in which particles spend a long time trapped by a surrounding cage and thereafter they escape from it. In fact, dynamically slowed systems begin to exhibit a plateau in the self-correlation function, reflecting the time spent in cages, and then a decay ($`\alpha `$-relaxation) reflecting escape from the cage and free configurational motion. This approach to the glass transition was formalized  with the idea that motions in dense fluids can be divided into intra- and inter-basin motions, where the basins refer to the multidimensional potential energy as a function of the particles coordinates. This separation of time scales, and thereby type of motion, is reasonable for dense systems, and may be used to justify a conceptual partitioning of the entropy into two parts, a configurational and a local contribution . These ideas can be re-expressed by considering the system, at a fixed average density, composed of central particles trying to escape their cages of neighbours, which are themselves fluctuating, and exchanging with their neighbours. Evidently, the $`\alpha `$-relaxation process corresponds, in systems very close to arrest, to the escape of the central particles from their cages after some time. In the absence of attractive forces, or for relatively long-ranged ones, the motions of the cage are restricted by the packing forces. This picture is relatively clear, at least in a phenomenological manner, for hard core spherical particles. For systems that possess a strong repulsion, and a long-ranged attraction, there are no new special features. Thus both “central” and “cage” particles remain in their mutual range of attraction while the structural rearrangement takes place. The resulting cage breaking can be viewed as being almost the same as that for hard spheres, but with changed zero of energy. Even if we acknowledge that there are effects due to attractions, they can still be considered to be weak perturbations of the picture arising from the hard core spherical particles.

Now let us turn from the scenario just described, typical of hard-sphere particles or particles where the attractions are long-ranged compared to the core size, to situations where attractions play a principal role. In this case, the freedom of the cage particles is considerably reduced. Indeed, particles must remain within a certain distance from each other where substantial attractive energies are still available, or they lose the advantages of being in the liquid-like structure. This being so, the cage around the central particle is much more rigid, since only smaller excursions from average positions of these cage particles are possible. In fact, for sufficiently narrow wells, the time that a particle spends inside the cage increases; a plateau regime results, and the system will eventually freeze. This “cage rigidity” was also the determining factor in hard core particles at very high density, leading to the typical colloidal “repulsive glass”. Here, though the mechanism and detailed laws will be different, we see that short-ranged attractions are able to cause cage rigidity, and formation of a solid, either a glass, or a crystal. Following these arguments we see that for particles with short-ranged attractions, the density and temperature window over which the liquid state is stable is greatly reduced since the configurational entropy is reduced. The system separates into a low-density fluid state where there is sufficiently free motion, and a dense state where there is in any case not sufficient configurational entropy to sustain a liquid, and the system freezes. This state will be a crystal or a dynamically arrested state. Based on the developing views of others, and what we shall present here, all these expectations are borne out from precise calculations, and increasingly in experiments.

From the discussion above, we may expect that when attractions are short-ranged the arrested glass is of different nature than the typical repulsive one, since it is favoured by both the energy (because the particles are not close enough to sample mainly repulsive energy), and the local entropy (given that the density is yet relatively low). We have earlier called this arrested state an attractive glass , and the crystal, previously found by a number of others , we name the attractive crystal. This distinguishes it from the typical face centred cubic (FCC) crystal formed by repulsive forces and the analogous repulsive glassy state.

Given the possibility that for short-ranged attractions, solids, both crystals and glasses, can be formed by these two distinct mechanisms leading to cage-rigidity, we may suppose that it is possible in principle for them to co-exist. For fixed (short) ranged potential one way of changing the balance of attraction and repulsion is to change the density, and we comment that such co-existence has been previously shown for crystals , and glasses , and is reproduced here.

In summation, we expect that, as the well width is narrowed, the liquid state becomes progressively less favoured, and is replaced by a conventional repulsive crystal, or its equivalent glass. An attractive crystal, and its equivalent repulsive glass, should also be present at higher densities, both types of solid co-existing at some typical densities where attractions and repulsions compete. From the preceding considerations we expect that these predictions should be general, irrespective of the detailed shape of the potential, and reflecting the typical range of the potential.

We have studied two potential energy models typically used to mimic colloidal interactions, the square well (SW) and the hard-core Yukawa potential, using a variety of techniques of condensed-matter theory. We have determined the phase diagrams calculating the liquid free energy using perturbation theory for the SW and the self-consistent Ornstein-Zernike approximation (SCOZA) for the Yukawa potential. The crystalline free energy has been calculated applying second order perturbation theory for both potentials. This technique has been applied to short-ranged potentials before , with quite remarkable success in reproducing the crystalline free energy in comparison to Monte Carlo experiments. However, it has been noted  that the same method is not so satisfactory in calculating the free energy of the liquid and gas states for very narrow attractive ranges, and this affects the accuracy of the phase diagram . To avoid this problem in the case of the Yukawa potential, we have modified the calculation so that the proven good features of the perturbation theory of the crystal are combined with the proven quality of the SCOZA for the liquid and gas free energies to produce a phase diagram that is of uniformly good quality. Currently there is no working SCOZA method for the SW case, though there is no fundamental barrier to developing one.

In our work the ideal Mode Coupling Theory (MCT)  of super-cooled liquids has been used in order to locate the glass transition curves. For colloids this method has also been found to describe many elements of the transition to the arrested state , though it does suffer from some limitations, a matter to which we return later. In the MCT calculations we have used the structure factors from Percus-Yevick approximation (PYA) for the case of SW and from SCOZA for the case of the Yukawa potential. These have proved to be quite accurate theories of the liquid and fluid states. The SCOZA method in particular has been compared to Monte Carlo simulations for a range of screening parameters (“well-widths”) of the Yukawa potential, and it has been shown  that the agreement for the phase diagrams is quantitative, at least for modest values of screening parameter.

By combining the results from the different methods described above, we are able to give, for some regimes, what we believe to be quite accurate phase diagrams. For the general case we believe that the results are at least qualitatively correct, and provide us a coherent picture of the connection between well-width and the arrangement of gas, liquid, crystal and glass phases. It is this overview of how the various phenomena fit together that is currently missing, and that should prove useful in the various applications alluded to above. We also point out that, from experiments, it has become clear that the interaction between globular proteins, in the range where they may crystallize, is characterized by short range attractions . In the case of proteins, a better comprehension of the phenomena would imply the study of anisotropic types of potential , due to hydrophilic-hydrophobic patching of the protein surface. The usual approach is, however, to use an effective isotropic interaction obtained by averaging over the anisotropy. Thus, as discussed later, most of our conclusions will be also important to understand the process of protein crystallization.

The paper is organized as follows. In the next section II we describe the approximate closures to the Ornstein-Zernike liquid integral equation that we used for the Yukawa and SW potentials. The methods employed to determine the equilibrium fluid and solid phases are described in sections III and IV. Section V is devoted to a brief sketch of MCT applied to attractive potentials, while in section VI the relevance of the spinodal curve for colloidal systems is discussed. The complete phase diagram, including the structural arrest lines is discussed in section VII, while section VIII is devoted to our conclusions.

## II Approximations to the OZ equation for the Yukawa and Square Well models

In this section we shall discuss the theory used in the investigation of the phase diagram. The Ornstein-Zernike (OZ) equation for the pair correlation function $`h(r)`$ is

$$h(r)=c(r)+\rho 𝑑𝐫^{^{}}c(|𝐫𝐫^{^{}}|)h(|𝐫^{^{}}|)$$
(1)
where $`g(r)=h(r)+1`$ is the radial distribution function and $`c(r)`$ the direct correlation function. Another important quantity is the static structure factor $`S_q`$, which is the equal time correlation function of the density variables in wave vector space,

$$S_q=\rho _q(t)\rho _q(t)/N$$
(2)
where the average $`\mathrm{}`$ is performed at equilibrium and the density variables are $`\rho _q(t)=_ie^{i𝐪𝐫_i(t)}`$, where the sum runs over all $`N`$ particles in the system. The Fourier transform of the correlation function $`h_q`$ is related to this quantity by the relation $`S_q=1+\rho h_q`$. The OZ relation in the wave vector space reads

$$S_q=\frac{1}{1\rho \widehat{c}_q}$$
(3)
$`\widehat{c}_q`$ being the Fourier transform of the direct correlation function.

As it stands Eq. (1) is not closed and some type of approximation is needed in order to solve it. In what follows we have chosen to calculate structural and thermodynamical properties using the PYA for the SW potential and the SCOZA  for the Yukawa potential. Both these model potentials possess some fundamental properties that make them good candidates for studying properties of attractive colloidal particles when the range of interaction is short. The use of SCOZA has been justified by the success of such approach in predicting simulation data.

### A The SCOZA for the Yukawa potential

The application of SCOZA to a hard-core Yukawa fluid provides a semi-analytic calculation of the thermodynamic properties of the fluid, liquid and gas states of the system . SCOZA has been applied to Yukawa systems with relatively large values of the range of the potential, with a satisfactory reproduction of the liquid-vapour binodal curves and a good description of the critical point region.

The hard-core Yukawa fluid is described by the following inter-particle potential,

$$v(r)=\{\begin{array}{cc}\mathrm{}r<\sigma & \\ \sigma ϵ\frac{\mathrm{e}^{b(r\sigma )}}{r}r\sigma & \end{array}$$
(4)
The parameter $`ϵ`$ defines the energy scale, while the parameter $`b`$, known as screening parameter, determines the range of the potential. The larger the $`b`$, the shorter is the range of the potential. In this paper we set $`\sigma =1`$ and $`ϵ=1`$, therefore the screening parameter is in units of the reciprocal of the hard-core diameter, and the temperature in units of the well depth.

SCOZA provides a closure relation for the OZ equation (1) by expressing the direct correlation function $`c(r)`$ in terms of the potential $`v(r)`$, as for other approximations like the mean spherical approximation (MSA), PYA, hypernetted chain (HNC) approximation etc., but introducing one or more state dependent parameters that can be adjusted to force the system to satisfy various exact thermodynamic relations of the system. In particular, the simplest SCOZA studied so far assumes

$$g(r)=0r\sigma $$
(5)
while $`c(r)`$ for $`r\sigma `$ is composed of two contributions, describing respectively the soft part of the potential and the hard core, as

$$c(r)=A\beta v(r)+K_{HS}\frac{\mathrm{e}^{b_{HS}(r\sigma )}}{r}r\sigma .$$
(6)
The Yukawa function in Eq. (6) takes into account the contribution to $`c(r)`$ arising from the hard-core repulsion. Thus, the two parameters $`K_{HS}`$ and $`b_{HS}`$ can be determined by setting $`v(r)=0`$ in Eq. (6) and requiring that both the compressibility and the virial route to thermodynamics lead to the Carnahan-Starling equation of state for a hard-sphere fluid . This amounts to describing the hard-sphere correlations via the Waisman parameterization . The soft part of the contribution in (6) is assumed to be proportional to $`v(r)`$ and hence it has the same range of the potential. The proportionality constant $`A`$ is calculated by imposing the condition that the compressibility and the energy routes yield the same result. This corresponds to the condition ,

$$\frac{}{\beta }\widehat{c}(q=0)=\frac{^2}{\rho ^2}\left(\frac{U^{ex}}{V}\right)_T$$
(7)
where $`U^{ex}`$ is the excess internal energy . Eq. 7 implies a partial differential equation (PDE) for $`A(\rho ,\beta )`$. The Yukawa potential (4) lends itself particularly well to the SCOZA scheme, because for this kind of interaction it is possible by means of Eqs. (5), (6) to establish an analytic relation between $`\widehat{c}(q=0)`$ and $`U^{ex}`$, which allows us to obtain straightforwardly a closed PDE from Eq.(7) by using $`U^{ex}`$ as the unknown quantity instead of $`A(\rho ,\beta )`$ . Once the internal energy has been obtained by numerical solution of this PDE, the Helmholtz free energy is calculated by integration with respect to $`\beta `$.

The reasons which led us to adopt SCOZA can be summarized as follows. The method uses a reasonable choice of functional relationship between $`c(r)`$ and $`v(r)`$ on the basis of numerous calculations obtained over a number of years. Furthermore, it has been shown  that the results obtained are in excellent agreement with computer simulations, at least for not too narrow ranges of the potential. In this sense it may be viewed as the best semi-analytical method to study the Yukawa potential. Finally, it is simple and convenient to obtain numerical solutions of Eq. (7), and this is helpful to make a survey of a problem in a large parameter window, rather than just in a small part of the parameter space. To give the reader some sense of the accuracy of SCOZA in some of the regimes we will discuss later, we present in Fig. 1 and 2 the phase diagrams for $`b=1.8`$ and $`b=6`$ calculated using SCOZA, compared with Gibbs Ensemble Monte Carlo (GEMC) studies . The data by Shukla are the largest system sizes with this potential yet studied by GEMC. In all cases, the agreement with SCOZA is excellent, and we may, in this regime of $`b`$-values and in calculating the phase diagrams, consider SCOZA to be equal to the best simulations. The method is superior to the other simple closures, including for example MSA.

Having said this, we note that there is little real information available about the detailed reliability of the closure relation where the potential becomes much narrower than for $`b=9`$, though it is reasonable to suppose that many properties are still satisfactory for somewhat larger $`b`$-values. Another comment we may make is that, when the range of the potential narrows, Eq. (6) may not be the optimal closure to ensure that $`c(r)`$ is accurate. Both these points should be regarded as words of caution, and as potential directions to develop the SCOZA method to better represent the structure in such problems. Indeed this will be the subject of future work.

### B The Percus-Yevick closure for the SW potential

We shall now briefly discuss the closure for a fluid of colloidal particles interacting via a square well potential,

$$v(r)=\{\begin{array}{cc}\mathrm{}r<\sigma & \\ ϵ\sigma <r<\sigma +\delta & \\ 0\sigma +\delta <r& \end{array}$$
(8)
where, in the present discussion, we set $`ϵ=\sigma =1`$ and we define the square well parameter $`\mathrm{\Delta }=\delta /(\sigma +\delta )`$, which parameterizes the attractive range of the potential. This model has been already the object of great interest in colloidal science . The state of the system is specified by three control parameters, the packing fraction $`\varphi =\pi \rho \sigma ^3/6`$ ( where $`\rho `$ is the number density, i.e. $`\rho =N/V`$), the temperature $`k_BT`$, and the square-well parameter $`\mathrm{\Delta }`$ of the attractive shell.

The PYA for $`c(r)`$ is $`g(r)=0`$ for $`r<\sigma `$ and

$$c(r)=g(r)\left[1\mathrm{e}^{\beta v(r)}\right]$$
(9)
outside the hard core . We solve the OZ equation in PYA using Baxter’s method of the Wiener-Hopf factorization . This correspond to rewriting Eq. (1) in terms of the real factor function $`Q(r)`$, defined for $`r>0`$. For $`0rR`$, $`R`$ being the range of the potential ($`\sigma +\delta `$ in the present case), one has

$$rc(r)=Q^{}(r)+2\pi \rho _r^R𝑑sQ^{}(s)Q(sr)$$
(10)
as well as, for $`r>0`$,

$$rh(r)=Q^{}(r)+2\pi \rho _0^R𝑑s(rs)h(|rs|)Q(s).$$
(11)
$`Q(r)`$ determines $`S_q`$ via its Fourier transform,

$`S_q^1`$ $`=`$ $`\widehat{Q}(q)\widehat{Q}(q)^{},`$ (12)
$`\widehat{Q}(q)`$ $`=`$ $`12\pi \rho {\displaystyle _0^{\mathrm{}}}𝑑r\mathrm{e}^{iqr}Q(r).`$ (13)
The resulting equations, obtained implementing the PYA in Eqs.(10) and (11), are then solved numerically to calculate the structure factor $`S_q`$ (see reference  for more details).

## III The Perturbation Theory Applied to Solid and Fluid phases

### A The Crystal Phase

In this section we shall discuss the method used to calculate the solid free energy. The perturbative approach has previously been used by Gast et al.  to construct the phase diagram of a colloidal solution with depletion interactions.

Further applications of the method appear in references . For example, in  the authors compared the results obtained for both phase diagrams and free energies of the solid and liquid phases for an Asakura-Oosawa potential and made comparisons to Monte Carlo calculations. They noted that the phase diagrams are overall reasonably satisfactory, but that the crystal free energies are excellent, in most regimes being essentially quantitative. It was concluded that the potential limitations in accuracy of the phase diagrams arise from the use of perturbation theory to the fluid and liquid phases, rather than the crystal. Also, for the crystal free energies there does not appear to be a significant loss of accuracy when the range of the potential becomes narrow. We have, therefore, applied perturbation theory to the crystal state of the Yukawa potential. Besides what we show here, it is opportune to note that there are in the literature other perturbative approaches, such as in .

The method is summarized as follows. We separate out explicitly the interaction potential as a hard-core contribution plus the attractive tail. The hard-core part of the potential is used as a reference for the perturbation, and the attractive tail is the perturbation itself . In other words, we decompose the potential as,

$$v(r)=v_0(r)+v_{att}(r)$$
(14)
where $`v_0(r)`$ is the hard-core repulsive potential and expand around the reference state $`v_0(r)`$. With this choice, the zeroth order term of the free-energy expansion coincides with the hard-sphere free energy. Once the perturbation expansion is carried out to second order we have for the Helmholtz free energy the following expression ,

$`{\displaystyle \frac{\beta F}{N}}={\displaystyle \frac{\beta F_0}{N}}+{\displaystyle \frac{\beta \rho }{2}}{\displaystyle v_{att}(r)g_0(r)𝑑𝐫}+{\displaystyle \frac{\beta F_2}{N}}`$ (15)
here $`\beta F_2/N`$ is the second-order perturbation term. $`F_0`$ and $`g_0`$ are respectively the Helmholtz free energy and the radial distribution for the reference hard-sphere system.

We focus our attention on the second order term in the expansion in Eq. (15). Indeed, its exact evaluation requires the calculation of higher order distribution functions , which are very hard to compute or to approximate reasonably. Barker and Henderson  proposed an approximation to $`F_2`$, based on the following observation. Since,

$$\frac{\beta F_2}{N}=\frac{1}{2}\beta (W_N^2_0W_N_0^2)$$
(16)
where $`W_N=_{i<j}^Nv_{att}(|𝐫_i𝐫_j|)`$, Barker and Henderson proposed to divide the space into concentric spherical shells and to calculate averaged properties using the number of particles in each shell. Following this route, they rewrote Eq.16 in terms of average numbers in the shells,

$$\frac{\beta F_2}{N}=\frac{1}{2}\beta \underset{ij}{}(N_iN_jN_iN_j)v^iv^j$$
(17)
where $`N_i`$ is the number of particles in the shell $`i`$ and $`v_i`$ is the perturbation energy, considered constant, within the shell. The first approximation consists of ignoring the correlations between shells, i.e.

$$N_iN_jN_iN_j=0$$
(18)
for $`ij`$. Moreover, inside a given shell, a second approximation is made,

$$N_i^2N_i^2N_ik_BT(\rho /P)$$
(19)
The two approximations (18) and (19) are equivalent to considering the volume of the shells to have the compressibility properties of a macroscopic portion of space (for more details see ). As a result, Barker and Henderson approximated the second order term in the expansion as,

$$\frac{\beta F_2}{N}=\frac{\beta \rho }{4}\left(\frac{\rho }{P}\right)_0v_{att}^2(r)g_0(r)𝑑𝐫$$
(20)
This approximation was found to be satisfactory in all calculations carried out so far. In our work the integrals in Eqs. (15) and (20) have been performed by a five-point integration rule, while for differentiation a central-difference scheme has been used. .

To carry out the calculation, we require the Helmholtz free energy and radial distribution function of the unperturbed hard-sphere system in the solid phase. It has been shown by computer simulation that a hard-sphere fluid shows a solid-fluid transition, for which the fluid phase alone exists up to a packing fraction $`\varphi =0.49`$ and the solid FCC phases exists for $`\varphi >0.55`$ . In between, there is a two phase coexistence of solid and fluid. These properties are well studied and the various information required in perturbation theory can be deduced from these studies. We note in passing that recently a renewed interest has been arisen in the equilibrium structure of a hard-sphere crystal. Indeed it has been believed for a long time that hard spheres crystallize with an FCC structure. Confocal microscopy observations, however, have rather found a random hexagonal phase which consists of a stack of FCC and hexagonal close packing (HCP) layers . Simulation seems to explain the phenomenon in terms of the small free-energy difference between FCC and HCP structures . In this paper we assume the crystal equilibrium structure to be FCC since this is believed to be the more stable .

To provide continuity with previous authors we make the choices described below. The state equation for a hard sphere FCC solid has been proposed by Hall  who derived a phenomenological expression based on computer simulation results, i.e.

$`Z_{HS}={\displaystyle \frac{P_{HS}V}{Nk_BT}}=`$ (21)
$`{\displaystyle \frac{1+\varphi +\varphi ^20.67825\varphi ^3\varphi ^40.5\varphi ^56.028\varphi ^6f(\varphi )}{13\varphi +3\varphi ^21.004305\varphi ^3}}`$ (22)
with $`f(\varphi )=\mathrm{exp}((\pi \sqrt{2}/6\varphi )[7.93.9(\pi \sqrt{2}/6\varphi )])`$. The compressibility can be derived by differentiating the compressibility factor $`Z_{HS}`$ as,

$`\left({\displaystyle \frac{\rho }{p}}\right)_0={\displaystyle \frac{\beta }{Z_{HS}+\varphi \left(Z_{HS}/\varphi \right)}}`$ (23)
In order to calculate the excess hard-sphere Helmholtz free energy $`F_0^{ex}`$ from the compressibility factor $`Z_{HS}`$, a thermodynamic integration in the packing fraction $`\varphi `$ can be performed, obtaining

$$\frac{\beta F_0^{ex}}{N}(\varphi )=\frac{\beta F_0^{ex}}{N}(\varphi ^{})+_\varphi ^{}^\varphi \left(Z_{HS}1\right)\frac{d\eta ^{}}{\eta ^{}}$$
(24)
Since the zero-density limit of a FCC crystal cannot be represented as easily as the one for a gas, alternative routes to perform the thermodynamic integration in (24) have to be devised . We have chosen to perform the integration starting from a packing fraction value of $`\varphi ^{}=0.544993`$, for which the value for the free energy has been calculated by computer simulation to be $`F_0^{ex}(\varphi ^{})N=5.91889`$ . We recall that the excess Helmholtz free energy is defined as the excess with respect to the ideal gas contribution .

For the radial distribution function $`g_0(r)`$ we have used an analytic formulation proposed by Kincaid and Weis that fits Monte Carlo simulation for a hard sphere FCC solid . This formulation is known to provide a good estimate of the hard sphere radial distribution function, at least in the range $`0.52\varphi 0.56518`$. Eqs. (15) and (20) can now be solved.

Once the Helmholtz free energy is evaluated, following the route we have just described, the Gibbs free energy and the pressure can be calculated as

$`\beta G`$ $`=`$ $`{\displaystyle \frac{\left(\rho \beta F\right)}{\rho }}`$ (25)
$`\beta P`$ $`=`$ $`{\displaystyle \frac{\rho \beta G}{N}}{\displaystyle \frac{\rho \beta F}{N}}`$ (26)
We have earlier noted that the perturbation theory for the crystal is highly accurate. Indeed, the second order perturbation term is useful, but it is interesting to note that the great bulk of the free energy correction for narrow well problems is captured by the first order term alone. To understand this, it is worth reflecting on the fact that the free energy of the crystal in the presence of short-ranged attractions is referred to the hard core crystal, and there is no question of the perturbation theory having to determine a priori any gross structural information. The corrections from attractions arise by virtue of the small changes in local vibrations that the particles make around their lattice positions, a portion of these motions involving the particles being within their mutual attractive range. In first order perturbation theory these contributions are treated as if the nature and distribution of the vibrations is unchanged, and the additional attractive energy contributions calculated essentially as an integral over the attractive potential multiplied by the zeroth-order hard core correlation function. The fact that second order contributions are typically small is suggestive. In a heuristic manner, we may argue that the intrinsic limitation on the extent and complexity of the phase space of the localized particles, imposed by their being in a crystalline state, means that even when attraction is incorporated, the changes in the nature and distributions of these vibrations is small. We note that this rationale is clearly inapplicable to the case of liquid, gas, and fluid states, where the addition of attractions to significantly affect the distribution of particles motions. Possibly this is the reason why perturbation theory works well for the crystal, indeed far outside its expected limitations, but is less successful for the other states. For completeness we note that for the crystal (Yukawa and SW) the errors, as estimated by the ratio of second-order to first-order terms are typically of order $`0.5\%`$. For the liquid they can be larger. We have however, also studied a square-well fluid using perturbation theory, and used these results along with those for the crystal to generate a phase diagram. The perturbation theory of the square-well fluid is therefore briefly discussed below.

### B Liquid phase for the SW fluid

In this section we shall discuss the method we adopted to calculate the thermodynamical properties of a fluid of colloidal particles interacting via a SW potential.

We chose the hard-sphere fluid as the reference system and treated the attractive part as the perturbation. The natural choice to describe the thermodynamics of a hard-sphere fluid is the Carnahan-Starling (CS) equation of state ,

$$\frac{\beta P}{\rho }=\frac{1+\eta +\eta ^2\eta ^3}{(1\eta )^3}$$
(27)
The CS equation provides an accurate account of the thermodynamic behaviour of the hard-sphere fluid for the entire region of the fluid phase. Its very simple analytical form makes it possible to obtain a closed expression for the Helmholtz free energy by integrating over density, as in Eq. (24). The zero-density limit of the free energy is the ideal gas value, so the thermodynamic integration starts from zero density. Thus, we obtain

$$\frac{\beta F^{ex}}{N}=\frac{\eta (43\eta )}{(1\eta )^2}$$
(28)
The compressibility is evaluated as in Eq. (23) by differentiation. For the radial distribution function we use a modification of the analytical PYA $`g_0(r)`$ for a hard-sphere fluid that was proposed by Verlet and Weis  to overcome certain limitations of this closure. Indeed, in PYA , the contact value $`g_0(r=\sigma )`$ of the radial distribution function underestimates the real value obtained by computer simulation and, also, the oscillations of the tail are slightly out of phase and too weakly damped. Verlet and Weis proposed,

$$g_0(r/\sigma ,\varphi )=g_0^{}(r/\sigma ^{},\varphi ^{})+\delta g_1(r)$$
(29)
Here, $`g(r)`$ is written as the sum of two terms. The first term corresponds to the solution for $`g_0(r)`$ within PYA, but evaluated at a smaller packing fraction value $`\varphi ^{}`$ and a smaller diameter $`\sigma ^{}`$, while the second term is a short-ranged correction $`\delta g_1(r)`$. The parameter $`\sigma ^{}`$ is then evaluated via a minimization of the difference between the simulation result for $`g_0(r)`$ and the PYA analytical result, between $`1.6\sigma `$ and $`3\sigma `$. This contribution improves the long-range behaviour of the PYA result. The addition of the short-range term $`\delta g_1(r)`$ improves the value at $`r=\sigma `$. Analytical forms for $`\varphi ^{}`$ and $`\delta g_1(r)`$ are given by the authors . The improved radial distribution is within $`1\%`$ of the computer simulation result in the whole range of packing fractions. With this result, the Helmholtz free energy for the hard-sphere reference system is calculated, and we may then proceed as in the previous section to calculate the Gibbs free energy $`G`$ and pressure $`P`$ using Eqs. (25) and (26).

## IV The Construction of Phase Diagram of the Yukawa Model by Hybrid Method

In the case of the Yukawa potential we have modified the approach of previous researchers somewhat in order to obtain the benefits of the best methods of condensed and liquid state theory. We have used SCOZA to calculate the liquid, gas and fluid phases free energy, but applied perturbation theory for the crystal free energy. To mark the difference with previous calculations, where phase equilibrium lines had been calculated by perturbation theory both for the crystal and for the fluid phase, we name our approach a hybrid method.

Phase boundaries between two phases (gas-liquid, fluid-solid etc.) are obtained by imposing the standard conditions,

$`\mu ^{(1)}`$ $`=`$ $`\mu ^{(2)}`$ (30)
$`P^{(1)}`$ $`=`$ $`P^{(2)}`$ (31)
where $`\mu `$ is the chemical potential, i.e. the Gibbs free energy per particle, $`\mu =G/N`$.

In those cases where SCOZA and perturbation theory are quantitatively validated the equilibrium phase diagram is highly accurate. As noted above, and shown in Figures 1 and 2, the SCOZA is well validated up to values of $`b=9`$. For values of $`b`$ less than or equal to $`6`$, results are indistinguishable in terms of phase equilibria from the best simulations that have been carried out . Similarly, as we shall discuss later, the perturbation theory rarely produces an error of more than $`0.5\%`$ in the free energy of the solid phase, although this analysis is based on certain assumptions about the perturbation series. Combining these observations, we believe that our phase diagrams are quantitatively accurate up to at least $`b=9`$. Beyond that, we make no particular claim, except that we expect that this hybrid method should still remain superior to the typical theoretical approximations that have been applied previously. Simulations have not been carried out beyond $`b=9`$.

## V Glass Transition and Mode Coupling Theory

### A Theory

The study of the glass transition in colloidal systems has been one of the most striking cases of verification of the current theories of super-cooled liquids. Early experimental studies  involved colloidal particles that are very closely represented by hard spheres where only excluded volume effects are important at high concentrations. Moreover, in contrast to simple atomic liquids, it is possible to avoid the crystalline phase beyond volume fractions of $`49\%`$ for sufficiently long periods of time to study the glassy-type dynamical processes, and ultimately the colloidal glass. The agreement between certain aspects of MCT and experiments on colloids is quite satisfactory  and the details of the time correlation functions are quite well reproduced. It is widely believed that, in deeply supercooled molecular liquids, the slow dynamics involves more complex dynamical processes than those described by MCT, and there the theory becomes of more qualitative applicability. Thus, the case of colloidal particles is of some practical interest in applying this type of theory.

In fact, even for colloids, small discrepancies appear in the comparison between experiments and MCT. The most important is the value of the critical volume fraction for the hard-sphere arrest transition, the experimental value being about $`58\%`$, while MCT predicts about $`52\%`$. This is of little importance where the dynamical laws at the hard-sphere transition are being compared between experiment and theory. Previous researchers have applied a shift to the transition volume fraction, and then fitted the laws in this region . Since the only current information on arrest driven by attractive interactions is that provided by MCT , there is as yet no accepted manner in which we can correct the MCT curves. This is somewhat inconvenient in the current context since for some parts of the parameter space the equilibrium phase-diagrams are quantitatively accurate, and it would be very satisfying to be able to superimpose, without correction, the relevant MCT arrest curves.

We now briefly review the nature of MCT, and discuss the type of information it yields. The MCT of super-cooled liquids describes the non-ergodicity transition by a nonlinear integro-differential system of equations for the normalized time correlation functions of density fluctuations $`\mathrm{\Phi }(q,t)`$. Apart from parameters entering from the microscopic motion, the only input to the MCT equations is the equilibrium wave-vector dependent structure factor of the system, $`S_q`$. The glass transition lines can be identified by studying the long time limit of the MCT equations, which determine the non-ergodicity parameter of the system $`f_q=lim_t\mathrm{}\mathrm{\Phi }(q,t)`$. An ergodic state is characterised by $`f_q=0`$. This value is always a solution of the MCT long-time limit equations . Thus, the glass transition appears as an ergodic to non-ergodic transition for the system, where $`f_q0`$ solutions arise. These points, thus, correspond to bifurcation singularities of the MCT equations, and, depending on the number of control parameters of the model, these can be of increasingly higher order, producing interesting features of the arrested states diagrams.

A good quality $`S_q`$ is an important input for a good description of the MCT arrest transition, as for the equilibrium phase diagram. In the earliest discussions of colloidal systems with short-ranged attractive interactions the Baxter interaction , a limiting case with an infinitely deep and zero-ranged SW attractive potential, was discussed by a number of authors . Subsequent studies indicated that the MCT equations are pathological for this interaction . The calculations were therefore extended to a SW potential both in the PYA and the MSA . Another solution of the MCT equations was obtained using the Yukawa potential and the MSA . Some common aspects emerged in these works. However, the SW model at first appeared to give a richer behaviour for the arrest transition curve. Thus, the system was shown to possess a glass transition curve in the parameter plane $`(\varphi ,k_BT)`$ the shape of this curve depending on the value of the SW parameter $`\mathrm{\Delta }`$ . For narrow well-widths, two branches of the glass curve have been identified. These have been interpreted respectively as a transition between a fluid phase and repulsion dominated glass (this is the typical repulsive glass) and a transition between a liquid and an attractive-interaction dominated glass (named the attractive glass). The two branches join and for approximately $`\mathrm{\Delta }4.11\%`$, a glass-glass coexistence between the two different types of glass appears. This coexistence line terminates in an end-point, beyond which the non-ergodicity parameters become the same for the two types of structures. The relevant singularity points, such as the end-point or the point where the glass-glass transition line reduces to a single point, are identified with higher order singularities of MCT equations and lead to unusual logarithmic dynamical relaxation laws . The mechanical properties of the system have been also studied  and they reinforce this interpretation. Earlier studies of the Yukawa potential did not locate this glass-glass phenomenon , but it was subsequently realized that the screening parameters which had been studied were not large enough, and further calculations seem to give clear indications that both the SW and Yukawa potentials give the same typical behaviour , implying that this does not crucially depend either on the potential shape or the approximation used for calculating the structure factor. Thus, it is now believed that this glass-glass scenario, and the attendant dynamical laws, is essentially a universal feature of the very short-ranged attractive potential.

We note in passing that the formation of two solid glass phases for very short-ranged potentials should not be too surprising. In fact we have earlier alluded to the fact that there are two crystalline phases in the phase diagram of such potentials. We may typically view glasses or arrested states as long-lived meta-stable states of the system that have not been able to equilibrate to the nearby crystal, and that are trapped in a restricted portion of phase space. In this sense we may expect each crystal to have associated to it a particular glass-type. Since one of the crystals in our phase diagram is “attractive energy dominated”, and the other “repulsive energy dominated”, it is hardly surprising that there should be two types of glass, dominated by the two regimes of interaction.

A further comment on the relevance of these glass curves is that we may see them as more than simply ‘transition curves’. If one reflects more deeply on the nature of the equilibrium phase diagram, and the thermodynamic states present in them, we recognize that at a deeper level they are reflecting the fact that for those particular parameters, the phase-space is dominated by a particular structure: a crystalline, liquid, or gas structure. The arrest curves carry analogous information. Thus, in the vicinity of the arrest curve, we may understand that most of phase space is becoming increasingly inaccessible, and breaking into smaller regions that are disconnected. That this may occur even for states that appear to have static structures typical of liquids is the distinguishing feature of glasses. Associated to this observation are those of dynamical slowing and tendency to arrest during any phase-separation through which such a glass curve passes. We shall argue in the conclusions that such phenomena are relevant in protein crystallization.

### B Experimental Studies

We now discuss some of the particular experimentally determined features that are associated with dynamical arrest driven by attractive interactions. For example, Verduin and Dhont  determined a curve of structurally arrested states in the phase-diagram of a system with short-ranged attractive depletion interactions. This locus, in some cases, intersects the binodal line and is referred to as transient gelation when observed in the spinodal region. Interestingly, it was these experimental authors who were the first to comment that MCT might also be applicable to cases where attractive interactions are important. As we noted above, subsequent results of such calculations have been most interesting . In addition a number of other experimental programs involving particles with depletion interactions have been published which offer many interesting insights for example in . In particular, Poon et al  have studied the arrest transition for systems where the range of the interactions is short, and their more recent work in this topic involves detailed connection to the theory described above . Other systems may have certain advantages over the depletion interaction system, but it is as yet too early to decide this issue . The results are typically quite promising, with some of these other systems also exhibiting some of the phenomena predicted by the theory. We may note in particular a more recent set of experiments that are intriguing in that they make detailed predictions for the correlation functions in a particular (re-entrant) part of the phase diagram. Thus, concentration time correlation functions have been observed in a polymer micelle system with a decay process much longer than the usual stretched exponential, and the results have been well-fitted to a logarithmic time-relaxation, , as predicted by the theory.

The development of experimental understanding, and deepening of the theory of systems with short ranged potentials is really just beginning, and many experimental programs have now been commenced or reoriented in efforts to make progress. However early information indicates that the MCT type theory may be able to describe main elements of the principal phenomena, at least for reasonably high volume fractions, where it is possible to separate aging from dynamical arrest in a reasonably clean manner.

## VI Are spinodals in Colloidal Systems meaningful?

Here we will take the liberty of raising a few issues in relation to spinodal curves that are calculated via SCOZA, or indeed many other typical liquid state theories. The reason that we make these comments is that such curves should have a particular status for these colloidal systems that are not relevant generally for molecular systems. We note first that the spinodal curve is determined from the condition that the curvature of the free-energy with respect to the relevant density variable becomes zero, and that this corresponds to the fluid phase becoming unstable as we lower the temperature. Between the binodal and the spinodal, the liquid and gas states in co-existence may be the global free-energy minima, but the fluid state remains meta-stable. Inside the region bounded by the spinodal, only the liquid and gas states in coexistence are stable. Now it is well known  that the free energy is a convex function, and it possesses only one minimum, and for some years now it has been understood that the spinodal curve determined from approximate theories (eg. mean-field theories) that consider two separate branches of the free energy and then connect them, has no real scientific basis. Careful Monte Carlo simulations carried out in systems of increasing size have lead to the conclusion that the spinodal curve shifts with the system size, merging with the binodal in the limit of infinite systems. There is some loosely defined kinetic phenomenon however, though even there it is not possible to define a spinodal curve, but a cross-over regime where the kinetic mechanisms begin to change from nucleation and growth to more collective phenomena. Interestingly enough, when the particle size becomes large these more sophisticated expectations are less relevant. Thus, it transpires that the relevant parameter in this story is the ratio of the particle diameter to the correlation length of the fluid. For very large particles such as high molecular weight polymers, colloidal and other particles, the microscopic length is so large that one has to be extremely close to the critical point to see fully developed fluctuations beyond the mean-field type ideas. Another consequence of this is that critical exponents in such systems as proteins and micelles have often been measured with mean-field values because experiments were not performed in the true critical regime. Similarly, the normal scepticism about the existence of a spinodal curve should be less relevant here, and we may expect the colloidal systems to exhibit quite reasonable spinodal behaviour. We have therefore included the spinodal curve in our phase diagrams.
In concluding this section on dynamically defined objects in the phase diagrams, it is commented that we have chosen to plot the MCT curves through the meta-stable regions between binodals and spinodals. Again, in colloidal systems, for the reasons given above, it is to be expected that such curves would have meaning, whereas they would not be meaningful in molecular systems.

## VII Phase-Equilibria and dynamical arrest lines

In this section we combine the results from the different techniques described above to exhibit the state of the system for given well depth.

### A Yukawa potential

In this section results for the Yukawa potential with the temperature in units of the well depth ($`k_BT/ϵ`$) plotted against the volume fraction of the system, $`\varphi `$ are presented. The values of the well width are determined via the screening parameter of the Yukawa $`b`$. The hard core radius is fixed at unity, so all quoted lengths are in units of the hard core diameter. We show results for values of $`b=5`$ (Fig. 3), $`b=6.05`$ (Fig. 4), $`b=7.5`$ (Fig. 5), $`b=30`$ (Fig. 6), $`b=60`$ (Fig. 7), and $`b=100`$ (Fig. 8). For comparison we also show (Fig. 9) a calculation for the square well system with $`\mathrm{\Delta }=0.03`$.

We begin with the largest well width, corresponding to $`b=5`$, see Fig. 3. Here the well width, considered say as the distance of half amplitude of a Yukawa, is comparable to the particle hard-core size. This is the typical situation that we are familiar with in elementary phase diagrams of atoms and molecules where van der Waals interactions predominate. Thus, we see the expected pattern of phase behaviour. Below the critical temperature, the gas-liquid phase equilibrium occupies the greater part of the low- and middle-range of densities, above the triple point. The crystal is favoured at higher density, and the liquid- and fluid-crystal boundary is nearly vertical, that is at fixed density, reflecting the substantial absence of any energy scale in the problem. The crystal is so tightly packed, and the attractions are so spread out across the system that it is only the repulsive part of the potential that is fundamental for crystallization. Indeed this is one of the important ideas in traditional liquid state theory; that attractions are not relevant to crystallization.

The asymptotic limits of solidification and melting boundaries of this coexistence at high temperatures reflect the hard-sphere limits respectively of $`\varphi 0.49`$ and $`\varphi 0.55`$, as expected from many simulations and theoretical observations . The triple point temperature is labelled $`T_p`$. Within the crystal phase there is also the boundary for the dynamically arrested state (the MCT transition line), again with high-temperature asymptote of $`\varphi 0.52`$(circles), the hard-sphere volume fraction for the glass transition predicted by MCT. We recall that MCT underestimates the glass transition packing fraction by about $`0.06`$. To call attention on this shift we report the true MCT curve (circles) as well as the MCT curve shifted by $`0.06`$ in packing fraction using $``$ as symbol. This boundary is also almost vertical, again reflecting the fact that, for wider wells, the arrest transition is driven essentially, at high enough temperature, by the repulsive part of the potential. Thus, no attractive glass is observed for this range of the potential. It has been shown theoretically , experimentally  and by simulation , that on decreasing the range of the attractive potential the fluid-fluid coexistence curve becomes meta-stable with respect to the fluid-solid one. This means that for short enough potential ranges the liquid-gas critical point is hidden in a phase separating region. This has two main consequences. The coexistence curve between the fluid and the solid is very broad implying, for low temperatures, coexistence between a very low density fluid and a high density solid (FCC). Also, the critical point becomes meta-stable with respect to gas-solid co-existence. We will return to some of the possible implications of all this later in the discussion. We have calculated the particular value of the screening parameter, $`b`$, at which the critical point becomes meta-stable with respect to the low density fluid-gas equilibrium, to be $`b^{}=6.05`$, and then presented the phase diagram at that value (see Fig. 4). We note that slightly different values have been previously reported in the literature. Thus, Hagen and Frenkel quote the values $`b^{}=7.4`$ using a Monte Carlo perturbation theory, and the value $`b^{}=6`$ based on Gibbs Ensemble Monte Carlo (GEMC) . Their GEMC value is very close to the value we have found. However a number of other values have also been quoted in the literature. Menderos and Navascues  used a density functional approach to determine $`b^{}=8.25`$, whilst Shukla  quotes a much higher value of $`b^{}=13`$. This latter value in particular is much higher than previous ones and our value, and considering that it is based on more extensive simulations than previous research this might be a matter of concern. The issues in relation to the accuracy of SCOZA have been addressed for modest values of $`b`$ in Figures 1 and 2, where the liquid-gas phase diagrams were shown for $`b=1.8`$ and $`b=6`$. However, if we examine Table 4 and Figure 4 of , we can readily compare the simulations to the precise predictions from SCOZA for selected values of the vapour and liquid densities for a range of temperatures. For $`b^{}=6`$, we find remarkable agreement between SCOZA and the Shukla gas-liquid equilibria. Also, whilst they are not quoted, we may estimate the critical temperature for differing screening parameters from these same simulations, and interpolation of others, and again conclude that the discrepancy between what we find for $`b^{}`$, and the Shukla result cannot arise because of a difference between the SCOZA and GEMC results for gas-liquid systems. The problem lies in the estimation of the crystal free-energy, or entropy, or estimation of the triple point, the other aspect of what we need to know to determine $`b^{}`$. In  the freezing transition location is determined by the so-called one-phase entropic condition . This condition implies that the freezing density is essentially constant, and almost unchanged for the range of $`b`$ values 1.8 to 10. In fact, in Figure 4, the phase diagram for $`b^{}=6.05`$, we do see rather significant deviation of the freezing density in the vicinity of the critical temperature from its high temperature limit. Indeed, comparing the other phase diagrams, Figures 5-8, we conclude that this variation is intrinsic to the whole short-ranged scenario, since it is the prelude to the splitting of the solid phase into two crystalline phases. The true underlying discrepancy arises because the entropy criterion implies that attractive forces are not important in the regime of crystallization currently under discussion, whereas the perturbation theory implies that they are highly significant. Both approaches are approximations, but it is possible to estimate the errors arising from the perturbation theory by considering the higher order terms.

Thus, for $`b^{}=6.05`$, (here $`k_BT_c^{}/ϵ=0.454`$ and $`\varphi _c^{}=0.230`$) we can propose to bound the errors in the perturbation theory of the crystal by reporting the ratios of the second- to first-order terms in the perturbation series. In the regime where the low-density fluid and crystal are in equilibrium we find that the ratio of the second- to first-order terms is approximately $`0.0044`$, whilst the first order-term has an absolute value of $`10.4`$ . If this ratio represents a true estimate of the errors, then the perturbation theory would appear to be quite satisfactory, the curvature of the freezing curve genuine, and the attendant flattening of the fluid side of the coexistence also quite accurate.

In general by estimating the impact of such errors in perturbation theory, we can estimate the shift of the fluid side of coexistence, and thereby estimate errors. In so doing we find the typical error in $`b^{}`$ will be less than $`1\%`$ from this source. If we accept this means of characterizing the error in perturbation theory, a matter clearly based on the series being well-behaved at successive orders, then we might conclude that the value of $`b^{}=6.05`$, in agreement with Frenkel’s original calculation, is a good approximation. If this were confirmed, then the calculations reported here for the overall phase diagram are probably amongst the most accurate for modest values of the range parameter, despite the fact that simulation is never used. This is not the primary motivation of our paper, but it would be an interesting way of approaching phase diagrams in future. Despite these optimistic estimates, further more careful evaluations are required by different methods to find a truly accurate value.

We should not imply that the value of $`b^{}`$ is of such crucial importance in the overall picture offered here. However, it does provide a useful check between different researchers and methods of approximating the phase behaviour in that its accurate estimation requires some satisfactory and simultaneous treatment of gas, liquid and solid phases.

Finally, we note that ten Wolde and Frenkel have made some interesting comments in relation to the kinetic processes that might be expected in this regime .

In Fig.5 the case $`b=7.5`$ has been plotted. As noted by a number of authors, the meta-stable gas-liquid phase equilibrium curve has now flattened considerably, and the low-density fluid-crystal coexistence (on the fluid side) has nearly the same slope. Increasing curvature of the crystal side of the fluid-crystal equilibrium curve is observed, arising from the increased influence of attractions on the crystal. The binodal line and the spinodal line in this case lie completely within the region of fluid-solid phase separation, the triple point has disappeared and the critical point is buried below the flat part of coexistence curve between the low-density fluid and the solid. As we mentioned earlier, this meta-stable behaviour would normally not be observable due to fluctuations, but for colloidal systems, globular proteins and other large particles, we may expect to observe such phenomena. Thus, on quenching such a system we might expect to see rather a rich pattern of behaviour, depending on the density that we quench at, and the depth of the quench. In particular, it is noted that we should be able to see a meta-stable gas and a liquid, the latter arresting into a glass at sufficiently low temperature, because the glass curve crosses the binodal, and spinodal at a finite temperature. We have earlier alluded to the idea that critical fluctuations can play an important role in the formation of crystals, for example protein crystals . The present screening parameter regime would exemplify this type of phenomenon since here we have a meta-stable liquid and gas that are critical (for large particles this should have some observable lifetime). The equilibrium phase diagram exhibits a fluid-crystal co-existence, so we have the possibility of crystal nucleation and growth phenomena in the presence of this meta-stable critical fluid, and it is this matter that ten Wolde and Frenkel have discussed . However, we also note the point that for this value of the screening parameter, as yet, the glass curve has not begun to severely interfere with the gas-solid equilibrium curve, and this is an additional advantage in the formation of crystals. This is in distinction to subsequent phase diagrams where the glass curve extends across much of the space.

It is in fact worth reflecting on the shifted, more realistic placement of the glass transition curve (stars). Thus, we see that at high temperature, as expected, from $`55\%`$ to $`58\%`$ volume fraction values, we have a crystalline state that is not interrupted by a glass transition. However, as the temperature is lowered, the increased importance of the attraction leads to the glass curve crossing the low density fluid-crystal coexistence region, and beneath this, the crystal may never form without the glass being an alternative long lived state. This comment is relevant also to the case of $`b=6.05`$, but for $`b=7.5`$, the curvature of the solidification curve has increased greatly, so this effect is emphasized. This dramatic interruption of the crystallization scenario will become more and more significant as the range of the potential narrows, and will be an important theme in our discussion.

Now we turn our investigation to the case of very narrow potential ranges. In Fig. 6 we present the case $`b=30`$. The two-phase coexistence of fluid and solid now occupies a much larger portion of the parameter space. The gas-liquid critical point is hidden well below the fluid-solid coexistence line and on the crystal-fluid phase boundary the effects of critical fluctuations will be much less. It is interesting to note that the crystal side of the coexistence curve exhibits a strong deviation towards higher densities. Thus the coexisting solid will be much more dense that in the previous cases, since the strong short ranged attraction is pulling the particles closer at low temperature. This phenomenon is the precursor of a solid-solid phase coexistence that in this case is still meta-stable and lies within the sublimation curve.

The glass transition curve is most interesting. Of course, for $`T\mathrm{}`$, the curve asymptotes to the hard core value. The attractive forces at low temperature now begin to strongly affect the curve so that it now turns sharply to the left, passing very close to the submerged critical point. Thus, we may tentatively assign $`b^{}=30.0`$ as that value of the screening parameter at which the submerged critical point becomes submerged by the glass, as well as by the crystal-gas curve. The inherent inaccuracies in the MCT estimate, alluded to before, may mean that the $`b^{}`$ value may not be very accurate. However the phenomenon is interesting. It means that the meta-stable critical fluid is now competing with a glass transition. It will transpire that the nucleation rate of these fluid-crystal equilibria are very low, possibly due to the high interfacial tension and therefore it is feasible that one may be able to approach the liquid gas equilibrium and its nearby glass transition without significant interference of the crystal. The type of slowing that would arise from a combination of critical slowing down and glassy slowing down has not been discussed in the literature previously. It would be an interesting problem. However, the phenomenon discussed here would significantly affect the possibility to form high-quality crystals, perhaps rendering it essentially impossible. Even though the glass may eventually decay in favour of the crystal, this will never lead to high quality crystals. Of course, an alternative view of this situation is that under these conditions it may be possible to make interesting materials that have critical fluctuations frozen into the glassy phase. Protein scientists wish to make good quality crystal, materials scientists often wish to make interesting materials. Our comments are applicable to both situations.

It is interesting also to note that in this case the glass curve passes close to the gas-liquid critical temperature,but then dips and intersects the binodal line below the critical density, passing through the spinodal region. Such a scenario has been found by Verduin and Dhont  in experiments on colloidal systems. We should stress that for low densities the glass curve itself may not be reliable, as was discussed in . The situation in relation to this point is, as yet, not settled.
We now discuss the case of extremely narrow wells. It is worth noting that the previous two values of the screening parameters represent the typical range of values accessed by those studying depletion-induced attraction between colloidal particles  or globular proteins . The next set of phase diagrams (Fig.7 and Fig.8) represent the limit of these types of interaction, and may correspond to cases such as grafted coatings on latex particles where we can access much narrower ranges of potentials.
In Fig. 7 and Fig. 8 we present the cases $`b=60`$ and $`b=100`$ respectively. For numerical reasons, it is difficult to extend fluid-solid coexistence curves in the low-density region of the phase diagram. This is evident from the truncation of the phase diagrams, and has no fundamental significance. As already noted in this paper, SCOZA has not been tested with simulations for such narrow ranges of the attractive potential. On the other hand, the main phenomena appear using other methods, so at least these are expected to be reliable. We return to this point at the end of this section.
In both figures we observe similar features. The spinodal line is now buried deep in the sublimation curve as is the critical point. A most striking feature of these phase diagrams is the presence of a solid-solid coexistence. This first-order phase transition was already present for less narrow ranges but there it was meta-stable (see for example the shoulder in the phase diagram for $`b=30`$). This phase boundary represents the coexistence of two crystals with the same lattice structure but different lattice spacing and consequently different density. It is terminated by a critical point of the solid-solid coexistence. The origin of this co-existence is interesting. The presence of short-ranged attractive interactions causes competition with the hard-core repulsive interaction. The fact that both have variations that occur on very short length scales means that the system may be forced to ‘choose’ between the attractive-dominated crystal and the repulsive-dominated crystal.
In the low-density crystal region, i.e. $`\varphi <0.65`$, the crystallization is dominated by entropic effect. In other words the the system chooses to optimize the entropy to form an FCC structure. Increasing the density the particles become closer and at some density they are forced to remain in the attractive shell of their nearest neighbours. When this happens, there is a decrease in energy which leads to an “attractive” crystal. It is the energy that stabilizes the phase. It is indeed clear that such a phenomenon can be present only if the range of the attractive potential is short enough. We may note that an isostructural phase transition has been already discussed theoretically for other kinds of potential characterized by a short-range potential  and indeed was also detected by simulation . It is almost certainly a genuine phenomenon. Here the isostructural phase transition is present for both $`b=60`$ and $`b=100`$. Decreasing the range moves the critical point of the transition to higher density, and indeed this is also in agreement with the behaviour in simulations . It is interesting to note that it is possible to find a triple point $`T_p`$ at which the two solids and the fluid coexist at the same temperature.
In Fig. 7 and Fig. 8 the glass lines have been also plotted. They both tend to the hard-sphere limit for high temperatures, and bend towards low densities with decreasing temperature, as we have seen in the previous cases. For short enough interaction ranges the glass transition line does not pass close to the critical point. We note that for $`b=100`$ the glass curve appears to break into two branches, with an apparent discontinuity at that point marked $`P_D`$ in the figure. The low-density branch is called attractive glass, whilst the right-hand branch is called repulsive glass. For very short ranges, in other types of attractive potentials, we have located a glass-glass transition, a transition between two different type of glasses originated either by repulsion or by attraction . For the Yukawa fluid such a phenomenon also appears to be present, although it has not yet been investigated in detail .
Both the results for the Yukawa potential considered here and those for the square-well fluid that we present next for comparison clearly show that the distinction between the attractive solid and repulsive solid becomes sharper as the range of the potential becomes narrower.
To present in a coherent way the role of the attractions on the crystal, glass and liquid-gas coexistence line, we show in Fig. 10 the dependence of $`T_c`$ on $`b`$ and both the glass transition temperature and the solid-fluid first-order transition temperature at the critical packing fraction $`\varphi _c`$. For completeness, the inset shows $`\varphi _c(b)`$.

### B Square-well potential

Finally we discuss a single example for the SW model. By so doing we wish to make the point that the main phenomena that have been discussed above are independent of the details of the shape of the potential, and are essentially universal. We note, however, that the SW phase diagram is not expected to be so quantitatively accurate as that of the Yukawa potential for reason discussed in section III.
The SW model was solved as discussed in section II B. We discuss a case where the range of the potential is very narrow, i.e. $`\mathrm{\Delta }=0.03`$, and the result is presented in Fig.9. It is clear that the situation is very similar to that for the Yukawa fluid, i.e. a solid-fluid phase coexistence extends from high temperatures (where again it reaches the correct hard sphere limit) expanding dramatically towards low and high densities for low enough temperatures. An isostructural solid-solid phase transition, with a critical point, is also present in this case. Indeed a limited correspondence between different potentials based on their general characteristics (effective core, range and energy scale) has been recently proposed , and it is possible that this idea may have more general applicability. We hope to return to this more general concept of corresponding states at a later point .
The glass line also has a very similar shape to the Yukawa fluid. As noted above, and in earlier publications, we have a glass-glass transition, terminating in a higher order glass singularity,the A3 transition . We believe that the presence of these two types of glass, the attractive and repulsive glass, is the disordered analogue of the presence of the two types of crystal discussed in some detail above. In this case we may draw also an analogy between the presence of the isostructural critical point and the presence of this MCT singularity $`A_3`$ end-point at which the two glasses become identical. Also as the well width gets larger, the crystal-crystal critical point vanishes. We believe that the glass-glass analogous of this is the $`A_4`$ point .
It is quite reasonable to suppose that for every crystal, there should be an analogous glass, and for every critical point of such a crystal-crystal equilibrium there should be such an MCT singularity. It will be interesting to explore this idea in future. It has the appeal of a potential general joint classification of equilibrium and glass transitions.

## VIII Conclusions

We shall use the conclusions section of this paper for two purposes. Thus we shall attempt to sum up the practical conclusions of our calculations, but, at the same time, try to considerably broaden the discussion to make contact with the main experimental situations where they might be useful. So far, we have focused the discussion quite strongly on the narrowly defined consequences of studying a short-ranged hard core attractive Yukawa potential, so we will begin by summing up that aspect of the discussion.
Firstly, from the technical point of view, we have achieved a certain success combining a good liquid-state method of calculation with perturbation theory. The resulting hybrid method takes free-energies from SCOZA for the gas and liquid states, and from perturbation theory for the crystal. We have indicated, in broad terms, how this overall strategy could be applied to phase diagrams in general, and how it may be qualified by checking of errors and relating these errors to shifts in the phase boundaries. Given the potential to exploit powerful methods of liquid-state theory, and the remarkable success of perturbation theory for the ordered state when one chooses the correct zeroth-order state, this may be a competitive manner in which to proceed for many problems in future.
On the other hand, the methods to determine dynamical arrest (eg. MCT) are not nearly so developed, despite their relative success in colloidal science. In particular, the absolute values of the density and temperature at which the arrest takes place is not correct. This should not be surprising. Equilibrium theory has had the benefit of many more years of development, and much more effort devoted to bring it to this level of achievement whereas for the dynamical arrest alternative routes or different approaches have been developed only recently . However, this aspect is quite inconvenient.
As we try to apply this theory to more realistic situations in colloids, materials and biology, we see increasingly this important motif of a competition between the equilibrium Boltzmann view of the matter, and the dynamically arrested aspects. So far, these two fields have developed somewhat in independent manners. However from the point of view of these practical topics in nature, there is no distinction, and it is often not the separate behaviours, but the interplay and competition between them that is primary to the scientific issue. This paper is one of the first attempts to make connections between these phenomena, but one can clearly see the limitations. For greater insight into these topics, we will have to address the possibility to both improve systematically the methods of studying arrest transitions, and also their consistency with equilibrium transitions. This must represent one of the important technical challenges in coming years.
From a broader perspective, we have shown that when the range of the potential becomes short in comparison to the core size, the subtle interplay between entropy and energy begins to change its character. The range of densities over which configurational entropy is relevant is much reduced, and one begins to lose the liquid state, in favour of crystals, or arrested glassy states. The reasons have been discussed at the beginning of the paper. In essence they amount to the fact that to retain the benefits of short-ranged attractions the particles must not depart too much from their typical inter-particle distance, or they are no longer in their mutual attractive well. This loss of freedom of motion, and restriction of favourable configurations, leads to a lowered configurational entropy. Another way of expressing the same idea is that the short-ranged potential leads to the loss of easy fluctuations that can open the cage of neighbouring particles that trap a central particle. The probability of finding such an ‘opening’ of the cage is much reduced, and the time during which a particle is localized by its neighbours increases, diverging at the arrest transition to form the ‘attractive glass’ that we have discussed. As the range of the potential narrows, the means of egress permitted to the particles is further limited, and the attractive glass becomes more favoured. This glass is therefore an effective competitor to the liquid and crystalline phases of the system, and this is reflected in the fact that the fluid phase is eventually erased by the glassy phase, and the critical point submerged underneath the curve of arrest transitions. This aspect should not be confused with the equilibrium phase diagram, although it is interesting that there are many parallels between the two.
As we discussed before, when the range of the potential becomes very short, the competition between entropy and energy is responsible for the formation of two distinct crystal phases. Coming from the fluid side, a crystalline phase dominated by repulsion is present. This state arises from the fact that the, at such densities, the entropy of a system made of particles free to move only in their own Weigner-Seitz cell of a FCC structure is larger then that of the metastable fluid . We named this phase a repulsive crystal by analogy with the glass. On increasing the density, this crystal becomes unstable and makes a transition to a smaller Weigner-Seitz cell, so that particles mainly stay within the attractive well. The structure of this phase, that we named attractive crystal, is a more compact FCC crystal that optimizes the free energy by means of reducing the potential energy. The result at this point is that the crystal structure may adopt two different states. The situation is analogous to the case of a gas-liquid phase separation where the gas (low density) phase optimizes the entropy and the liquid (high density) has a much lower energy.
The implications of all of this are profound for practical situations. For reasons given earlier, many systems such as colloids, globular proteins, “buckyballs”, nano-particles, pre-ceramic particulates, and others have this property of a short-ranged attraction. In all these cases “precipitation”, “gellation”, “glassification”, or solidification are frequently the commonly observed outcome. In cases where we consciously seek to make such a state this is satisfactory, and it remains only to adjust the potential to have sufficiently short range to obtain the required properties of the solid. However, for cases such as globular proteins, and nano-scale or meso-scale ordered materials with prescribed optical properties the situation is quite different. Here we seek to make a crystal. In fact, reviewing the phase diagrams in Figs. 5-8, we can see why the crystal is hard to access. If we work to the right side of the glass curve (point marked $`A`$ in Fig. 6), universally we may expect to fall into the glass state; there is essentially no choice. Since the glass curve moves to low density, this is a serious restriction. However, we may choose to work within the two-phase low-density fluid-crystal co-existence regime, but to the left of the glass curve (point marked $`B`$ in Fig. 6). The outcome is then a question of kinetic control, and will not be completely settled by diagrams such as we are drawing. However we can make some educated comments. Thus, if we work in the two-phase regime of gas-solid ($`B`$), we may nucleate and grow crystals. Whether the proximity to a meta-stable critical point is advantageous or not, as discussed by Hagen and Frenkel  is not our primary concern here, though this is an interesting proposition. The broader point is that by nucleating to the left of the glass curve one may enter the crystalline region (thus form a macroscopic crystal) by a route not described by the “adiabatic” description here, and thereby avoid some of the complications of the glass. This renders the formation of crystal at least feasible, although where the glass curve runs through the two-phase region, it will remain difficult to form truly large high quality crystals. From this region where it is possible to crystallize, one should also exclude the two phase region, whether it is meta-stable or not, since the partial (micro-)phase separation, crystallization and glassification, all competing dynamically is unlikely to produce a good crystal also. This leaves only that region bounded to the right by the glass curve, to the bottom by phase separation, and to the top by the gas side of the gas-crystal phase-coexistence as a likely candidate for forming good crystals. This is interesting. It leads us to suppose that for a fixed short range of the potential there is a “practical crystallization region” in the temperature-density plane, irrespective of the specific features of the equilibrium phase diagram. However, more importantly, there is also a limited regime of interaction ranges where such a slot is significantly large enough to be accessed experimentally. We believe that much of the discussion that has taken place in the last few years in the literature in relation to protein crystallization  is almost certainly heading in the correct direction. Thus, Mushol and Rosenberg  show in their Figure 12 the typical situation for a globular protein phase diagram. They exhibit a phase diagram that has (a) two metastable liquid phases in equilibrium, a more and less dilute phase of protein (gas-liquid in our language), (b) a ’gellation curve’ that we associate with the glass curves in our work, and (c) a fluid-crystal coexistence regime. They name the ‘good’ regime for crystallization zone I, and the others II and III. We have also essentially partitioned our phase-diagram into the same types of zones, and concluded that this gas-crystal region described above by the glass, binodal and gas side of the coexistence boundary (zone I in their language) would be the most favourable for formation of crystals. For the sake of comparison in one of our figures (Fig.6) that has a range typical of globular proteins, we have marked regions I-III in analogy with ref..
We would argue that the present work, with its phase diagrams, and discussion of control parameters is the quantitative expression of these ideas that have arisen in the protein crystal literature. This is potentially encouraging, since it opens the possibility to make more quantitative study of these systems.
However, perhaps the most promising directions involve the study of the current model, and underlying ideas in more depth to see what independent kinetic routes exist to the formation of crystals. We have seen how these phase diagrams are indicative of the kinetic behaviour, but we believe that there will be much more significant insights as to how to form high quality crystals in these regimes if a deeper understanding of kinetics is acquired in future. In particular the glass analogy seems promising as a means to characterize the more confined phase space experienced by these systems. Aging, and kinetic phenomena in general, is an arena that is growing in importance  and may offer significant advances. The traditional viewpoint of activated processes, and simple kinetic processes is without doubt incomplete in the limit where we approach the rather confined phase-spaces characterized by approach to a glass transition. The realization that we are in a ‘glassy’ scenario may well assist in the development of new theories of kinetics of crystallization more appropriate for such questions.
In any case, one can hardly doubt the high degree of practical significance that kinetic phenomena associated with short-ranged potential systems will have in the coming few years. Given that we discuss a model potential that is only slightly different from those long considered in liquid-state theory, we must be intrigued by the novelty in supposedly simple situations.

###### Acknowledgements.

We acknowledge helpful comments by J.Bergenholtz and M.Fuchs who first pointed out to us that the Yukawa model would also possess an $`A_3`$ and $`A_4`$ scenario, but at much large screening parameter then their previous studies. Also, the most helpful advice and comments of R. Evans and D. Frenkel are gratefully acknowledged.
The work in Rome is supported by PRIN-2000-MURST and PRA-HOP-INFM, and the work both in Rome and in Dublin is supported by COST P1.
The work at Stony Brook was supported by the Division of Chemical Sciences, Office of Basic Engineering Sciences, Office of Energy Research, U.S. Department of Energy.

## Figures