# Stability of a hard-sphere binary quasicrystal

## §1. INTRODUCTION

It is well known that quasicrystalline phases are obtained from alloys of at least two metallic elements, mostly by processes of rapid solidification of the melt. However, the question of how quasicrystals form is far to be understood from the theoretical viewpoint (Steinhardt and Di Vincenzo 1991, Janot 1994). In this respect, an important step towards this goal could be given by a systematical study of numerical simulations at which a quasicrystalline phase is formed from a binary fluid. Unfortunately, this program is difficult to carry out since simulations of such systems are very expensive due to the large number of atoms which are necessary to display the quasiperiodic features. Recently, however, it was reported a first effort of a molecular dynamics simulation of quasicrystal solidification from the melt of a diatomic three-dimensional system (Roth, Schilling and Trebin 1995). In fact, it was shown that a binary Lennard-Jones fluid, at adequate cooling rates, indeed solidifies into a phase of icosahedral long-range order close to a quasicrystalline structure known as the truncated icosahedral binary model (TIBM). This is assumed to be a simplified version of the Henley-Elser model for the icosahedral structure of (Al,Zn)<sub>49</sub>Mg<sub>32</sub> (Henley and Elser 1986). A previous investigation by means of relaxation simulations had shown that the TIBM was at least strongly metastable for T=0, a feature which was mainly interpreted due to the fact that such structure comprises a very dense packing of spheres of two sizes (Roth, Schilling and Trebin 1990).

In the present paper we shall investigate the stability of a TIBM composed by hard-spheres by means of a density-functional theory (Lutsko and Baus 1990), which was previously applied to a one-component hard-sphere (HS) quasicrystal (Cataldo and Tejero 1995). HS solids may serve as suitable reference systems for mean-field- or perturbation-theory treatments of more realistic potentials (Daanoun, Tejero and Baus 1994, Achrayah and Baus 1998). Recently, Denton and Hafner (1997a,b) have studied a model of monatomic quasicrystal interacting via effective metallic pair potentials, by means of a combination of classical density-functional theory and thermodynamic perturbation theory, which takes as its reference system the above mentioned one-component HS quasicrystal.

Closely related to freezing and thermodynamic stability properties of HS solids is the problem of packing. This is a well studied subject for identical hard spheres, but less is known about packing properties of mixtures of spheres or disks of two different sizes, see e. g., Likos and Henley (1993).

The possibility of quasicrystalline metastable states for one- and two-component HS solids has been recently explored within density-functional theory in the range of diameter ratios above 0.85 (M<sup>c</sup> Carley and Ashcroft 1994). Then it was found that such phases have higher free energies than the crystalline or liquid states. We shall see that a similar conclusion stands for the TIBM for diameter ratios above 0.83. However, this trend is partially reversed for lower diameter ratios, since the quasicrystal becomes stable with respect to crystal but remains metastable with respect to fluid.

In the following Section we shall describe the TIBM sphere packing, showing that a small shift in the original position of large spheres makes the maximum packing fraction to become competitive with respect to HS crystal structures of similar composition. In Section 3 we shall give a brief account of the density-functional approach to the calculation of quasicrystal free energy. Finally, our results and conclusions will be presented in Sections 4 and 5, respectively.

## §2. TIBM SPHERE PACKINGS

Let us consider a decoration of the 3D Penrose tiling at which all vertices and middle edges of both rombohedral cells are decorated by identical spheres. Thus, if we take rombohedral edges as unity, the maximum diameter turns out to be 0.5. On the other hand, each prolate rombohedron can accomodate two larger spheres located symmetrically at both sides of the center of its long diagonal (see fig. 1). If the location of such spheres divides the diagonal into three parts of ratios $`\tau ^2:\tau ^3:\tau ^2`$, with $`\tau =(\sqrt{5}+1)/2`$, the maximum diameter (given by the separation of the centers of both spheres) is 0.563, yielding a diameter ratio 0.888 and a maximum packing fraction of 0.580. Here it is convenient to define a notation for the HS parameters namely, $`\alpha =\sigma _1/\sigma _2`$, $`x_2=N_2/N`$, and $`\eta =\frac{\pi }{6}[(1x_2)\alpha ^3+x_2]\frac{N}{V}\sigma _2^3`$, denote respectively, the diameter ratio, the concentration of large spheres, and the packing fraction.

Now, taking into account that in the Henley-Elser model the above ideal position of the larger spheres was chosen essentially for simplicity (Henley and Elser 1986), we may allow a greater separation between both spheres in order to increase the maximum diameter. In fact, by this procedure the maximum packing fraction can be improved to $`\eta _{max}=0.648`$, for a diameter ratio $`\alpha =0.798`$ when large spheres of adjacent rombohedra touch each other (in addition to the contact between both large spheres of each rombohedron).

In order to measure the relevance of the above packing fractions, we may consider a binary crystal of similar composition namely, a fcc lattice plus a decoration of small spheres on faces and large spheres on vertices, which has a concentration $`x_2=0.25`$ close to the TIBM value (0.236). In addition, we may compare to a disordered-fcc structure in which both kind of spheres are distributed irregularly over the sites of a fcc lattice. This disordered crystal allows us to select the exact composition of the TIBM.

In table 1 we give the packing fractions of the above structures for three selected values of the diameter ratio. We note that the diameter ratio 0.798 yields the better quasicrystal performance with a packing fraction slightly higher than the crystal one (the packing fraction of the disordered crystal is only relevant at the higher diameter ratio).

## §3. CALCULATION OF THE HELMHOLTZ FREE ENERGY

The above calculation of packing fractions may be regarded as an “order-zero” of a stability analysis of the competing HS structures. In order to get a much more definite analysis, in what follows we shall resort to a density-functional approach, the so-called generalized effective liquid approximation (Lutsko and Baus 1990), which has proven to provide a very accurate description of one-component HS solids. In particular, the application of this formalism to a monatomic HS quasicrystal of high packing fraction (Henley 1986), was recently reported (Cataldo and Tejero 1995), showing that it is thermodynamically unstable (i. e., it does not have a coexistence with the fluid), in agreement with previous studies (M<sup>c</sup> Carley and Ashcroft 1994). Thus, our present treatment will be based upon a straightforward generalization to mixtures of the above approach. In fact, following Xu and Baus (1992), we write the Helmholtz free energy per particle of a HS binary solid as the following functional of both partial local number densities $`\rho _i(𝐫)`$ ($`i=1,2`$),

$$f[\rho _1,\rho _2]=\frac{k_BT}{N}\underset{i=1}{\overset{2}{}}𝑑𝐫\rho _i(𝐫)[\mathrm{ln}(\mathrm{\Lambda }_i^3\rho _i(𝐫))1]+f_{ex}[\rho _1,\rho _2],$$
(1)   
where $`k_B`$ denotes Boltzmann’s constant and $`\mathrm{\Lambda }_i`$ the thermal de Broglie wavelenght of species $`i`$. The excess free energy per particle of the solid $`f_{ex}`$ is approximated by the one of an effective HS binary fluid of the same concentration, $`f_{ex}/k_BT\psi (\widehat{\eta })`$, whose packing fraction $`\widehat{\eta }`$ is obtained from the solution of a system of two coupled differential equations (Tejero and Cuesta 1993),

$$\widehat{\eta }^{}(\lambda )=\frac{z(\lambda )\psi (\widehat{\eta }(\lambda ))}{\lambda \psi ^{}(\widehat{\eta }(\lambda ))}$$
(2)   
$$z^{}(\lambda )=\mathrm{\Phi }(\widehat{\eta }(\lambda ))$$
(3)   
as $`\widehat{\eta }=\widehat{\eta }(\lambda =1)`$. In the above equations, the prime denotes function derivative, the initial conditions are $`\widehat{\eta }(0)=z(0)=0`$, and

$$\mathrm{\Phi }(\widehat{\eta }(\lambda ))=\frac{1}{N}𝑑𝐫𝑑𝐫^{}\underset{i=1}{\overset{2}{}}\underset{j=1}{\overset{2}{}}\rho _i(𝐫)\rho _j(𝐫^{})C_{ij}^0(|𝐫𝐫^{}|,\widehat{\eta }(\lambda )),$$
(4)   
where $`C_{ij}^0`$ denotes the Percus-Yevick approximate direct correlation function of a binary HS fluid mixture (Lebowitz 1964). The HS density profiles are parametrized in terms of Gaussians centered at the sites $`\{𝐑^{(i)}\}`$ occupied by species $`i`$:

$$\rho _i(𝐫)=\left(\frac{\alpha _i}{\pi }\right)^{3/2}\underset{𝐑^{(i)}}{}\text{e}^{\alpha _i(𝐫𝐑^{(i)})^2}i=1,2.$$
(5)   
Replacing (5) in (4), and performing the angular integrations, we obtain<sup>1</sup><sup>1</sup>1In case of a disordered structure, instead of (5) we must consider the alternative expression, $`\rho _i(𝐫)=x_i\left(\frac{\alpha _i}{\pi }\right)^{3/2}_𝐑\mathrm{exp}[\alpha _i(𝐫𝐑)^2](i=1,2)`$, where the summation runs over all lattice sites. Accordingly, the expression (6) has to be modified by the replacements, $`_{𝐑^{(i)}}^{N_i}x_i_𝐑^N`$.

$$\mathrm{\Phi }(\widehat{\eta }(\lambda ))=\frac{1}{N}\underset{i=1}{\overset{2}{}}\underset{j=1}{\overset{2}{}}\underset{𝐑^{(i)}}{\overset{N_i}{}}\underset{𝐑^{(j)}}{\overset{N_j}{}}_0^{\mathrm{}}𝑑rrC_{ij}^0(r,\widehat{\eta }(\lambda ))S(r,\alpha _{ij},r_{ij}),$$
(6)   
being

$$S(r,\alpha _{ij},r_{ij})=\left[\frac{\alpha _{ij}}{2\pi r_{ij}^2}\right]^{1/2}\left[\mathrm{exp}(\alpha _{ij}(rr_{ij})^2/2)\mathrm{exp}(\alpha _{ij}(r+r_{ij})^2/2)\right],$$
(7)   
where $`\alpha _{ij}=2\alpha _i\alpha _j/(\alpha _i+\alpha _j)`$ and $`r_{ij}=|𝐑^{(i)}𝐑^{(j)}|`$. The calculation of (6) deserves some discussion. First, we note that the integral in $`r`$ can be analytically computed and, in the case of periodic structures, the double lattice sum reduces to a simple one running over spherical shells of sites centered around an arbitrary site. On the other hand, the lack of periodicity of the quasicrystal does not allow this simplification, making necessary to resort to another method. In a previous work (Cataldo and Tejero 1995), it was proposed a method of evaluation of quasilattice sums that can be applied to the present case. We refer the reader for details to such work, but in a brief summary the method consists in a generalization of the above mentioned shell summations of periodic lattices to a situation in which each site produces a different “shell” sum due to variations in the surroundings, as is the case of a quasilattice. So, the double sum consisting in the average of such “shell” sums over an infinite number of sites is approximated by an average over the sites within a finite size quasicrystal “sample”. Here it is important to notice that in the case of a crystal, this procedure yields the correct value only if the “sample” possesses the exact composition of the infinite crystal (note that it is subjected to variations due to boundaries). So, we have adopted this exact composition prescription for the quasicrystal “samples”.

The Gaussian width parameters $`\{\alpha _1,\alpha _2\}`$ of (5) are determined so as to minimize the total free energy of (1). Fig. 2 shows the effect of quasilattice sums on the value of quasicrystal free energy; we see that it suffices to work with “samples” of $`N_2>1000`$ large spheres in order to reduce convergence errors to less than 0.1$`\%`$. In fact, the computation time needed for evaluating each ($`\alpha _1,\alpha _2`$) free energy value is around 11’ CPU in a workstation Alpha DEC 3000 for a “sample” of 4000 large spheres. Fortunately, the location of free energy minima in the $`\alpha _1,\alpha _2`$ plane is rather insensitive to the “sample” size, so we could work with small “samples” ($`N_2100`$) in locating minima and, at a second stage, a minimum value of free energy was computed for a large “sample” ($`N_24000`$) within a small grid of typically 20 ($`\alpha _1,\alpha _2`$) points.

An alternative treatment of the problem of quasilattice sums consist in replacing the quasicrystal by a suitable crystalline approximant (Krajčí and Hafner 1992, Goldman and Kelton 1993). Such structures are crystals with large unit cells where the arrangements of sites closely approximate the local structure of the quasicrystal. Instead of an aperiodic tiling of space by prolate and oblate rombohedra (Penrose tiling), approximants are based upon tilings by the same rombohedral units which are periodically repeated in the form of large cubic unit cells. Once such tiling is defined, the decoration of both kind of rombohedra with large and small spheres can be done in the same way as in the quasicrystal. Lattice sums over the approximant are then computed by taking any of its cubic unit cells as a “sample” in the above explained sense. It is clear that approximants having larger unit cells are expected to yield free energy values closer to the quasicrystal ones. However, for HS potentials only the most local structure in the neighbouring of each particle should be of importance, and then one expects that the smallest approximants would be adequate to obtain acceptable results. We have found that this is indeed the case; one of the smallest approximants is the so-called ‘1/1’ cubic structure (Elser and Henley 1985), which has 20 prolate and 12 oblate rombohedra per unit cell. It is slightly richer in large spheres (+0.84%) than the quasicrystal and thus has a 0.03% (0.11%) higher packing fraction for $`\alpha `$=0.888 ($`\alpha `$=0.798). Regarding its free energy, it turns out to be always slightly lower than the quasicrystal one, namely, for $`\alpha `$=0.888 the free energy differences vary from 0.13% (nearly within convergence errors of the quasicrystal values) at low packing fractions, up to 0.21% at high packing fractions, while for $`\alpha `$=0.798, the values increase from 0.21% up to 0.41%. From all these figures, it becomes clear that such differences in free energy are most probably due to the slightly higher packing fraction of the 1/1 approximant.

## §4. RESULTS

In fig. 3 we represent the free energy of the different phases for $`\alpha =0.888`$, note that the quasicrystal turns out to be always metastable with respect to the remaining phases. On the other hand, both crystals show an exchange of relative stability with respect to the fluid<sup>2</sup><sup>2</sup>2The fluid phase free energy was computed from the BMCSL equation (Boublik 1970, Mansoori, Carnahan, Starling and Leland 1971)., which is consistent with a freezing into the disordered solid as expected (M<sup>c</sup> Carley and Ashcroft 1994, Xu and Baus 1992).

Another feature of solid phases concerns to the degree of localization of the spheres at the corresponding sites. This is measured by the Gaussian widths which are shown in table 2. Note the particular behavior of small spheres in the quasicrystal, they present an appreciably lower localization than the large or small spheres of the remaining cases.

An additional test on quasicrystal stability may be performed by taking the distance $`D`$ between the centers of both large spheres of each rombohedron as a variational parameter. In fact, in fig. 4 we find free energy minima for $`D`$ values which are approximately equidistant from the original $`D`$ in the Henley-Elser model and the one that allows to maximize the packing fraction (both distances are indicated by dotted vertical lines). It is interesting to note that the increase of free energy at both sides of minima is accompanied by a loss of particle (small and large spheres) localization. Moreover, for each packing fraction there exists a critical maximum separation (star points in fig. 4) beyond which the free energy minimum in $`\alpha _1,\alpha _2`$ vanishes, yielding thus an end of mechanical stability. In order to understand such behavior, we first note from fig. 4 that the critical $`D`$ increases with $`\eta `$, and then for a fixed $`D`$, say 0.66, the quasicrystal becomes mechanically unstable below certain density, $`\eta 0.53`$ in this case. In other words, the bounds of mechanical stability in fig. 4 denoted as star points, were to be expected as arising from the familiar fact that any solid becomes unstable below certain density.

Let us now repeat our analysis for the lower diameter ratio $`\alpha =0.798`$. Fig. 5 represents the free energy of quasicrystal (for $`D=0.627`$), crystal, and fluid. We see that the improved quasicrystal packing fraction has contributed to stabilize such phase with respect to the crystal, but it remains unstable as compared to the fluid. Note that the crystal free energy is very close to the quasicrystal one at $`\eta =0.53`$, but unfortunately it was impossible to perform a reliable calculation of the crystal free energy for $`\eta <0.53`$ since the $`\alpha _1,\alpha _2`$ minimum at those low densities began to suffer an overlap with a spurious minimum arising from the flaw of Percus-Yevick approximate direct correlation function of the effective fluid at solid-like densities<sup>3</sup><sup>3</sup>3Such spurious minima were discussed by Lutsko and Baus (1990) for one-component HS crystals.. On the other hand, we have found that the crystal minimum in $`\alpha _1,\alpha _2`$ vanishes for $`\eta >0.59`$, this unexpected behavior seems to be caused by some kind of high-density flaw of our density-functional theory<sup>4</sup><sup>4</sup>4Recently it has been shown that the density-functional theory employed by M<sup>c</sup> Carley and Ashcroft (1994) is unable to predict the existence of certain high density, mechanically stable, HS solids (Tejero 1997,1998).. Nonetheless, the crystal density range shown in fig. 5 is enough to perform our comparison of relative stabilities. The Gaussian widths of both solid phases are shown in table 3, again we observe that the small spheres in the quasicrystal are less localized than the large ones, but both spheres are more localized than the corresponding values of table 2. Regarding the crystal, note the anomalous behavior of the last column; the loss of localization of large spheres at increasing densities finally yields a vanishing minimum in the $`\alpha _2`$ coordinate. Finally in fig. 6 we test again the quasicrystal stability against variations of the distance $`D`$ between large spheres. We observe features similar to those discussed from fig. 4 above, despite the shorter range of $`D`$ due to a longer large sphere diameter.

Concluding our results, we note that the exchange of relative stabilities between crystal and quasicrystal (cf. figs. 3 and 5) occurs at a diameter ratio 0.83 (see table 1) when the free energies of both solids differ within the order of quasicrystal convergence errors, along the whole range of packing fraction.

## §5. SUMMARY AND CONCLUSION

We have analyzed the stability of a HS binary quasicrystal by means of a density-functional method. We have first focused upon a TIBM structure with a diameter ratio 0.888 that maximizes the packing fraction. It was found that this HS structure turns out to be thermodynamically unstable i. e., it has a higher free energy with respect to crystals and fluid of similar composition, being a disordered crystal the solid phase of lowest free energy. Next we have studied the quasicrystal free energy dependence upon the distance between nearest large spheres, finding a minimum at a slightly more distance than the original TIBM one. In addition, we have found that a small increase of such distance allows to improve the maximum packing fraction to a competitive value with respect to a crystal of similar composition. In fact, it was found that for a diameter ratio around 0.8, the quasicrystal is stable with respect to such crystal, but it remains unstable as compared to the fluid.

In conclusion, as was pointed out by M<sup>c</sup> Carley and Ashcroft (1994), it appears that crystalline states are preferred for diameter ratios above 0.85 but, as was shown in the present paper, a small decrease of diameter ratio to 0.83 leads to a favored quasicrystalline state which, however, does not reach a coexistence with the fluid. Such thermodynamical unstability is most probably due to the assumed HS interparticle potential, in contrast to the Lennard-Jones one of the simulation performed by Roth et al. (1995).

## ACKNOWLEDGMENTS

The author is indebted to C. L. Henley and C. F. Tejero for helpful comments and to D. M. Jezek for useful discussions.

REFERENCES

ACHRAYAH, R., and BAUS, M., 1998, Phys. Rev. E, 57, 4361.

BOUBLIK, T., 1970, J. Chem. Phys., 53, 471.

CATALDO, H. M., and TEJERO, C. F., 1995, Phys. Rev. B, 52, 13269.

DAANOUN, A., TEJERO, C. F., and BAUS, M., 1994, Phys. Rev. E, 50, 2913.

DENTON, A. R., and HAFNER, J., 1997a, Europhys. Lett., 38, 189; 1997b, Phys. Rev. B, 56, 2469.

ELSER, V., and HENLEY, C. L., 1985, Phys. Rev. Lett., 55, 2883.

GOLDMAN, A. I., and KELTON, K. F., 1993, Rev. Mod. Phys., 65, 213.

HENLEY, C. L., 1986, Phys. Rev. B, 34, 797.

HENLEY, C. L., and ELSER, V., 1986, Phil. Mag. B, 53, L59.

JANOT, C., 1994, Quasicrystals (Oxford: Clarendon Press).

KRAJČÍ, M., and HAFNER, J., 1992, Phys. Rev. B, 46, 10669.

LEBOWITZ, J. L., 1964, Phys. Rev., 133, A895.

LIKOS, C. N., and HENLEY, C. L., 1993, Phil. Mag. B, 68, 85.

LUTSKO, J. F., and BAUS, M., 1990, Phys. Rev. A, 41, 6647.

MANSOORI, G. A., CARNAHAN, N. F., STARLING, K. E., and LELAND, T. W., 1971, J. Chem. Phys., 54, 1523.

M<sup>c</sup>CARLEY, J. S., and ASHCROFT, N. W., 1994, Phys. Rev. B, 49, 15600.

ROTH, J., SCHILLING, R., and TREBIN, H. R., 1990, Phys. Rev. B, 41, 2735; 1995, Ibid., 51, 15833.

STEINHARDT, P. J., and DI VINCENZO, D. P. (editors), 1991, Quasicrystals: The State of the Art (Singapore: World Scientific).

TEJERO, C. F., 1997, Phys. Rev. E, 55, 3720; 1998, Ibid., 58, 5171.

TEJERO, C. F., and CUESTA, J. A., 1993, Phys. Rev. E, 47, 490.

XU, H., and BAUS, M., 1992, J. Phys.: Cond. Matt., 4, L663.