# Freezing of polydisperse hard spheres

## Abstract

The fluid - crystal equilibria of polydisperse mixtures of hard spheres have been studied by computer simulation of the solid phase and using an accurate equation of state for the fluid. A new scheme has been developed to evaluate the composition of crystalline phases in equilibrium with a given polydisperse fluid. Some common assumptions in theoretical approaches and their results are discussed on the light of the simulation results. Finally, no evidence of the existence of a terminal polydispersity in the fluid phase is found for polydisperse hard spheres, the disagreement of this finding with previous molecular simulation results is explained in terms of the inherent limitations of some ways of modeling the chemical potential as a function of the particle size.

PACS numbers: 05.70.Fh, 64.70.Do, 82.70.Dd 

The phase behavior of polydisperse mixtures of hard spheres (PHS) has received some attention in recent years. Different theoretical approaches  and simulation methods  have been used to gain knowledge about the transition from a polydisperse fluid phase to crystal phase(s). The theoretical approaches  use to involve drastic approximations regarding the composition of the phases, and the results are often presented in form of stability diagrams , both facts should be carefully taken into account when interpreting the results. Bolhuis and Kofke  have studied the fluid - solid equilibria by using molecular simulation methods, finding a ”terminal polydispersity” in the fluid phase that they interpreted as the maximum polydispersity of a fluid which can originate a freezing transition, and related such a result with some experimental data. The origin of the ”terminal polydispersity” in Ref  will be addressed in this work.

Let $`P(\sigma )`$ be a given probability distribution function of particle diameters (PDFD). The distribution can be characterized by its moments, $`m_k=<\sigma ^k>/\sigma _0^k`$, where $`\sigma _0`$ is a reference diameter.

The thermodynamics of PHS fluids is very accurately described by the generalization of Salacuse and Stell  to the polydisperse case of the equation of state (EOS) due to Boublik and Mansoori, Carnahan, Starling and Leland (BMCSL) . In such an equation the pressure, $`p`$, can be written as: $`\beta p=\beta p(m_1,m_2,m_3,\eta )`$ where $`\eta `$ is the packing fraction: $`\eta =\pi Nm_3\sigma _0^3/(6V)`$. $`N`$ is the number of particles, $`V`$ is the volume. and $`\beta 1/(k_BT)`$, with $`k_B`$ being the Boltzmann’s constant and T the absolute temperature.

The excess chemical potential, $`\mu _{ex}`$, in the fluid phase takes the form: $`\beta \mu _{ex}\left(\sigma \right)=_{k=0}^3c_k\sigma ^k`$, where the coefficients $`c_k`$ depend on $`m_1`$, $`m_2`$, $`m_3`$ and either $`\eta `$ or $`\beta p\sigma _0^3`$.

The goal of the present work is to evaluate the fluid - solid equilibrium for a given PDFD in the fluid phase. This point of view is the main difference with the calculations of ref , however the statistical mechanics underlying both procedures is basically the same.

In order to study polydisperse systems is convenient to make use of the semigrand (SG) ensemble, where the pressure, the total number of particles, $`N`$, and the chemical potential differences between the different species and a reference one are fixed. For hard body interactions the basic thermodynamic differential relation reads,

$$d\left[N\beta \mu _0\right]=Vd\left(\beta p\right)+\beta \mu _0dN\underset{i0}{}N_id\left(\beta \mu _{i0}\right)$$
(1)   
where $`\mu _0`$ is the chemical potential of the reference species. The sum is done over the other components, $`N_i`$ is the number of particles of species $`i`$ and $`\mu _{i0}\mu _i\mu _0`$. Later, a continuous distribution of sizes will be used however the discrete description is kept, for the shake of clarity in the equations:

An imposed chemical potential distribution (ICPD) is used to perform the calculations, such a distribution should produce the required PDFD in the fluid phase. As a difference with the procedure in Ref  here the ICPD will depend on the pressure. We have taken advantage of the accuracy of BMCSL EOS. Such an equation let us to link the fluid phase composition with the chemical potential distribution at a given pressure.

Let $`P_0(\sigma )`$ be the expected PDFD, for instance, a Gaussian distribution centered at $`\sigma _0`$ and with standard deviation $`\sigma _0\lambda `$. We can use as input the values of $`\beta \mu \left(\sigma \right)`$ given by:

$$\beta \mu \left(\sigma \right)=\frac{\left(\sigma \sigma _0\right)^2}{2\sigma _0^2\lambda ^2}+\beta \mu _{ex}^{BMCSL}(\sigma ,\beta p,\lambda )$$
(2)   
The actual PDFD of the fluid, $`P(\sigma )`$ will be practically identical to $`P_0`$ due to the accuracy of the BMCSL EOS: The two contributions to the chemical potential can be grouped, by using an unique set of coefficients $`\{a_k\}`$:

$$\beta \mu \left(\sigma \right)=\underset{k=0}{\overset{3}{}}a_k\left(\frac{\sigma }{\sigma _0}\right)^k$$
(3)   
The coefficients $`a_k`$ will be functions of $`\lambda `$ and $`\beta p`$. The strategy to evaluate the fluid - solid equilibrium under the conditions stated above lies on Gibbs-Duhem (or Clausius-Clapeyron) integration schemes . The procedure is sketched as follows: Starting for a given point ($`\beta p_0,\lambda _0`$) in which both phases are in equilibrium, we obtain a trajectory on the $`(\beta p,\lambda )`$ plane that keep equilibrium conditions fulfilled. This can be done because the chemical potential differences can be written as functions of $`\beta p`$ and $`\lambda `$ through the coefficients $`a_k`$. Therefore, considering $`N`$ fixed:

$$d\left[N\beta \mu _0\right]=\left[V\underset{i0}{}N_i\left(\frac{\left(\beta \mu _{i0}\right)}{\beta p}\right)_\lambda \right]d\left(\beta p\right)+$$

$$\left[\underset{i0}{}N_i\left(\frac{\left(\beta \mu _{i0}\right)}{\lambda }\right)_{\beta p}\right]d\lambda $$
(4)   
In the limit of a continuous distribution of sizes we can write a Clausius Clapeyron analogue equation for the coexistence $`(\beta p,\lambda )`$ line:

$$\left(\frac{d\beta p}{d\lambda }\right)_{coex}=\frac{_{k=1}^3\left(a_k/\lambda \right)_{\beta p}\mathrm{\Delta }m_k}{\mathrm{\Delta }v_{k=0}^3\left(a_k/\left(\beta p\right)\right)_\lambda \mathrm{\Delta }m_k}$$
(5)   
where $`\mathrm{\Delta }`$ represents the difference between the values of the corresponding property in the two phases and $`vV/N`$. The values of the derivatives of $`a_k`$ with respect to $`\beta p`$ and $`\lambda `$ can be evaluated numerically.

The starting point in the Clausius Clapeyron integration (CCI) was the monodisperse hard sphere system ($`\lambda =0`$) where the equilibrium pressure is known  to be $`\beta p\sigma ^311.71`$. A second order predictor corrector has been used to advance in the integration. The integration step was $`\delta \lambda =0.0025`$, the initial slope was found to be zero. The fluid properties have been directly extracted from the BMCSL EOS. A number of tests for several points on the ($`\beta p,\lambda `$) trajectory were performed by carrying out SG Monte Carlo (SGMC) simulations on the fluid phase using $`N=256`$ and, within numerical accuracy, no differences between simulation and theoretical results were found. The solid phase was considered to be in a face centered cubic (FCC) ordering and its properties were evaluated by SGMC simulation.

Details of the simulation procedure will be published elsewhere . It suffices to say that three kind of moves were performed, i) translation of the spheres (following the standard procedures), ii) changes of a particle diameter, by choosing the new diameter with probability proportional to $`\mathrm{exp}[\beta \mu (\sigma )]`$ with $`\sigma [0,\sigma _{max}]`$, $`\sigma _{max}`$ depends on the hard sphere interactions and iii) Changes of volume, where we found convenient to scale simultaneously the size of the particles to enhance convergence on the sampling.

CCI were performed by systems with with $`N=108`$ and $`N=256`$. No significant finite size effects were found regarding the main conclusions of the work. Some results for $`N=256`$ are presented in figures. In Fig 1 the results of the pressure as a function of the polydispersity in the two phases in equilibrium. In figure 2 we plot the packing fraction of the coexisting phases as a function of the polydispersity of the fluid phase. In figure 3 we show how the average diameter of the crystal phase increases with the polydispersity of the fluid phase.

In some theoretical work  some conjectures have been made regarding the possibility of finding a fluid in equilibrium with two or more solid phases. This result, based on diagrams of phase stability, is not consistent with the phase rule, except in a number of singular points. The origin of such results lies in the theoretical approximation that the fluid composition (polydispersity) is equal to the overall composition in the solid phases(s). As can be seen in figures 1 and 3 such condition is not fulfilled except for very low polydispersities. However, one should not neglect the possibility of finding two (or more) crystalline phases which could enter into competition to become the solid phase in equilibrium with the fluid. The phase diagram of the binary mixtures of hard spheres  can be used as an example; In general for an equimolar fluid mixture $`x=1/2`$, the solid phase in equilibrium with the fluid is richer in the large component. As the size difference increases an eutectic point appears in the phase diagrams , it could happen that for given size difference the eutectic composition could become $`x_s=1/2`$, in that case we could have for larger size differences, that the stable solid phase could become composed mainly by the small spheres. In order to check such a possibility for the polydisperse system we have made some control simulations starting from FCC phases composed with particles smaller that $`\sigma _0`$ at pressures closed to the ones obtained in CCI. Those systems evolved either to produce the same solid appearing in the CCI or to the melting of the sample. These results seem to discard the change of stable phase in the procedure of increasing $`\lambda `$ (at least in the range studied in this work). It is clear, however, that such a competition between solid phases could appear as the freezing proceeds, the change of the composition of the fluid phase will alleviate the phase rule restrictions, and it is quite likely to happen for high pressures where the fluid could even disappear as an equilibrium phase. Other possibilities have not been considered here, for instance, the formation of crystal phases with a bimodal distribution of sizes, which could be accommodated, for instance, in a body centered cubic lattice.

Here we will discuss briefly terminal polydispersity in the fluid phase which appears in the results of Ref. . In that work the ICPD has the form:

$$\beta \mu _{i0}\left(\sigma \right)=\frac{\left(\sigma \sigma _0\right)^2}{2\sigma _0\lambda ^2}$$
(6)   
With this function a CCI scheme was performed from the monodisperse limit ($`\lambda =0`$), and it was found that $`m_10`$ when $`\lambda `$ increases and that the values of the reduced polydispersity, $`s_2=\sqrt{m_2/m_1^21}`$, in the fluid at equilibrium with a solid phase have to be less than a certain ”terminal polydispersity”, which the authors identified with some experimental results of crystallization of colloidal mixtures. In our simulations no such a terminal polydispersity appeared, however the results are roughly similar to those of Ref  until that point. This apparent anomaly can be explained in terms of the form of ICPD. In our scheme, the excess contribution to the chemical potential produces a positive value of the coefficient $`a_3`$ of the ICPD, favoring large diameters to compensate the effect of hard sphere repulsions, however the form of the ICPD given by Eq (6) produces a limit in the maximum packing fraction of a fluid phase with a certain polydispersity well below close packing. It is possible to estimate such a limit by computer simulation  or using the BMCSL EOS. In fact the terminal points in the diagrams of  correspond just to the crossing between the fluid branch of the phase diagram and the line of such a maximum packing as a function of the polydispersity. As pointed out in Ref  the end of the coexistence curve is conditioned by the ICPD, however such an end does not seem to correspond to any relevant physical situation.

The existence of a terminal polydispersity in the fluid phase has been also treated theoretically (See for instance Ref. ), however the results are strongly influenced by the restrictions to the size fractionation.

From the inspection of Fig. 1 one could think that the stability range of a polydisperse crystal with respect to the fluid can lie in a range between two pressures. This is not the case, we must emphasize that two points with the same value of the reduced polydispersity on the crystal ”branch” do not correspond to the same PDFD (See Fig. 2 and Fig.3). The one corresponding to the higher pressure is associated with a distribution with greater value of $`m_1`$ and a more negative value of the skewness  (i.e. the distribution has an asymmetric tail extending out towards small values of $`\sigma `$).

The authors acknowledge the financial support of DGICYT/Spain (Grant No. PB95-072-C02-02)