# Vacancy-stabilized crystalline order in hard cubes

^{a}Soft Condensed Matter, Debye Institute for NanoMaterials Science, Utrecht University, 3584 CC Utrecht, The Netherlands;^{b}University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, United Kingdom; and^{c}Heinrich-Heine Universität Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany

See allHide authors and affiliations

Edited by Athanassios Z. Panagiotopoulos, Princeton University, Princeton, NJ, and accepted by the Editorial Board August 22, 2012 (received for review July 23, 2012)

## Abstract

We examine the effect of vacancies on the phase behavior and structure of systems consisting of hard cubes using event-driven molecular dynamics and Monte Carlo simulations. We find a first-order phase transition between a fluid and a simple cubic crystal phase that is stabilized by a surprisingly large number of vacancies, reaching a net vacancy concentration of approximately 6.4% near bulk coexistence. Remarkably, we find that vacancies increase the positional order in the system. Finally, we show that the vacancies are delocalized and therefore hard to detect.

The free energy of crystal phases is generally minimized by a finite fraction of point defects like vacancies and interstitials. However, the equilibrium number of such defects in most colloidal and atomic/molecular crystals with a single constituent is extremely low. Exemplarily, for the face-centered cubic crystal of hard spheres, one of the few colloidal systems where the vacancy and interstitial fractions have been calculated, the fraction of vacancies and interstitials is on the order of 10^{-4} and 10^{-8}, respectively, near coexistence (1). As such, neither the vacancies nor the interstitials strongly affect the phase behavior, and so most studies of crystals ignore the effect of these defects. Nevertheless, vacancies and interstitials have a significant effect on the dynamics in an otherwise perfect crystal, as the main mechanism for particle diffusion is hopping of particles from filled to empty sites or between interstitial sites.

In this paper, we examine a system of hard cubes where, as we will demonstrate, one cannot ignore the presence of vacancies. Arguably, a cube is one of the simplest nonspherical shapes and the archetype of a space-filling polyhedron. Surprisingly, despite the simplicity of this system, we find that the stable ordered phase is strongly affected by the presence of vacancies, so much that vacancies actually increase the positional order and change the melting point. Remarkably, the fraction of vacancies in this system is more than two orders of magnitude higher than that of hard spheres or any other known typical, experimentally realizable, single-component atomic or colloidal system, reaching 6.4% near coexistence. Additionally, while purely hard (not rounded) colloidal cubes are yet to be realized, colloidal cubes with various interactions are now a reality (2⇓⇓⇓⇓–7), and it is likely that hard cubes will be realized in the future.

Here, we use Monte Carlo (MC) and event-driven molecular dynamics (EDMD) (8) simulations to examine in detail the effect of vacancies on the equilibrium phase behavior of hard cubes. The model we study consists of *N* perfectly sharp hard cubes with edge length σ in a volume *V*. Aside from hard-core interactions that prevent configurations of overlapping cubes, the particles do not interact. In both types of simulation (MC and EDMD), overlaps are detected using an algorithm based on the separating axis theorem (e.g., ref. 9). More information on the EDMD implementation for cubes is given in *SI Text*.

## Results

### Spontaneous Vacancy Formation.

The equation of state for a vacancy-free system of hard cubes has been the subject of a number of studies (10⇓–12) and was most recently examined by Agarwal and Escobedo (10). It clearly shows a single, first-order phase transition between a fluid and an ordered phase. The authors of ref. 10 identified the ordered phase at coexistence to be a liquid crystal mesophase (i.e., a cubatic phase), which is characterized by the presence of long-range orientational order along three perpendicular axes, but a lack of long-range positional order. However, the authors noted that finite-size effects made it difficult to determine the extent of the positional order in their system, and based this identification on their observation of finite diffusion in the ordered phase. At high densities, both refs. 10 and 12 agree that the ordered phase is a simple cubic crystal.

There is no fundamental reason why diffusion cannot occur in a crystal, hence this is an insufficient criterion for distinguishing between a cubatic phase and a crystal. To study more closely the range of the positional order along the ordered branch, we reexamined the intermediate density region (near coexistence) using highly efficient EDMD (8) simulations allowing us to access system sizes more than an order of magnitude larger than the ones considered in ref. 10. We simulated systems of *N* = 40^{3} = 64,000 particles starting from a simple cubic crystal lattice. Coexistence between the fluid and ordered phase is observed directly for overall packing fractions 0.455 ≤ η = *N*σ^{3}/*V* ≤ 0.480. Snapshots of typical configurations are shown in *SI Text*.

Looking in detail at the EDMD simulations for η between 0.52 and 0.56 (i.e., in the region where ref. 10 reported the cubatic phase), we noticed that in many simulations the crystal lattice spontaneously transformed in one of two distinct ways (see *SI Text*). In most cases, the simple cubic crystal did not maintain its original orientation in the box; instead, it rotated, introducing defects and frustrations to the crystal lattice. In others, the system spontaneously added extra layers (i.e., extra lattice sites). As long as the crystal lattice remains aligned with the simulation box, the number of lattice sites can easily be measured from the number of peaks in the three-dimensional density profile of the cubes. Fig. 1 shows a two-dimensional projection of such a density profile from simulations with *N* = 40^{3} particles. From this plot, we can determine that the system has *N*_{L} = 41^{3} lattice sites. Thus, the system spontaneously incorporated a large number of excess lattice sites into the crystal. The resulting crystal has a net vacancy concentration, of approximately α = (*N*_{L} - *N*)/*N*_{L} = 8%. It should be noted that because the volume and the number of particles in the system are fixed, it is generally not possible to reach an equilibrium concentration of defects in the system. However, the formation of extra layers and lattice distortions both significantly increase the number of lattice sites in the crystal, and suggest that the thermodynamically stable phase in this regime might be a vacancy-rich crystal structure. We note here that the systems sizes examined by ref. 10 would not have allowed for extra layers to form. Hence, we believe that the rotated, defective crystals we see are what the authors of ref. 10 identified as cubatic.

To examine the effect of vacancies on the crystal structure, we performed additional EDMD simulations on systems with various net vacancy concentrations α = (*N*_{L} - *N*)/*N*_{L} for a number of packing fractions η = *N*σ^{3}/*V*. The average global positional order in the system was measured using the positional order parameter averaged over all particles: [1]where **K** is a reciprocal lattice vector of the crystal under consideration and **r**_{j} is the position of particle *j*. In all our plots we have chosen **K** to correspond to a single lattice vector (i.e., with and *a* the lattice spacing). To set the net vacancy concentration, we placed *N* = (1 - α)*N*_{L} particles randomly on a simple cubic lattice with *N*_{L} = 40^{3} = 64,000 lattice sites, and then rescaled the volume (and thus the lattice spacing) of the box to the desired packing fraction η. Hence, the resulting system has a lattice spacing that depends on the chosen packing fraction η and net vacancy concentration α. These simulations show that there is a maximum in the global positional order as a function of net vacancy concentration α for varying packing fractions η (Fig. 2). An increase in order due to an increase in number of vacancies is unexpected because typically the presence of defects (such as vacancies) decreases the order. This observation suggests that adding defects reduces frustration in the crystal, potentially stabilizing a defect-rich crystal.

### Defect Realization.

The “net vacancy concentration” α, defined above, is simply the excess of lattice sites *N*_{L} compared to the number of particles *N*, divided by *N*_{L} (i.e., the fraction of lattice sites that does not have a particle associated with it). In a typical system (for instance, hard spheres), a vacancy is localized to a single lattice site and one can determine the number of vacancies by counting the number of empty lattice sites. In a hard-sphere crystal, a particle next to an empty lattice site is kept in place by its other neighbors. This is not the case for hard cubes, and, as a result, the way vacancies manifest in this system is very atypical. In hard cubes, entirely empty lattice sites are rarely seen even in the vacancy-rich crystals near coexistence. Instead, a defect manifests itself as a finite-length chain of particles along one of the three principal axes in the crystal, in which the particles have a slightly larger interparticle spacing than the average, as shown in Fig. 3. Hence, if a vacancy extends over four lattice sites, as is the case for one of the vacancies highlighted in Fig. 3, then the vacancy is realized by three particles sharing four lattice sites in a one-dimensional chain. Additionally, while a two-dimensional layer in a typical snapshot, such as in Fig. 3, shows regions of disorder, it should be noted that even at high vacancy concentrations, the crystal still shows a well-defined lattice spacing on average, which can be easily determined from the position of the peaks in the scattering function *S*(**k**) or the three-dimensional pair correlation function *g*(**r**) (Fig. 4).

It should be noted that the net vacancy concentration includes both vacancies as well as interstitials in the sense that each interstitial cancels a vacancy. However, the number of vacancies is higher than the number of interstitials resulting in the large positive net vacancy concentrations we find in this system. Similar to a vacancy, interstitials are also not localized in this system and occur by *n* particles sharing *n* - 1 lattice sites.

### Phase Diagram of Hard Cubes.

So far, we have established a relationship between the order in the system and the net vacancy concentration. However, there is no way to determine the equilibrium net vacancy concentration and the phase boundaries from the EDMD simulations. A proper determination of the equilibrium concentration of vacancies as well as the phase diagram requires free-energy calculations. The free energy per particle (*f* = β*F*/*N*) of the solid with vacancies is given by: [2]where the first term is the translational free energy of a normal Einstein crystal (13), the second term is the rotational free energy of the crystal (14), and the third term is the combinatorial entropy associated with placing *N* particles on *N*_{L} lattice sites: [3]Detailed descriptions on how these terms were calculated using MC simulations are given in *Materials and Methods*.

We determined the free energy as a function of net vacancy concentration α, at a fixed packing fraction η = 0.52. As shown in the inset in Fig. 5*B*, the minimum in the free energy occurs for a high concentration of vacancies. Specifically, we find that at this density the number of particles is 4% lower than the number of lattice sites. While calculating the free energy as a function of α, we also observed that the free energy, excluding the combinatorial contribution, was almost linear. This is shown in the inset of Fig. 5*B*. Although we are uncertain to the origin of this linearity, it seems to suggest that the vacancies are only weakly interacting.

For each value of α, the free energy as a function of density was obtained by combining the reference free energies shown in Fig. 5*B* with a separate equation of state measured for that value of α. By minimizing the resulting free energy with respect to α, we find that the number of excess lattice sites decreases as a function of the density, as expected (Fig. 5*C*).

Using a common tangent construction (see *SI Text*) in combination with the determined free energies, we mapped out the phase diagram that is shown in Fig. 5*A*. We find coexistence between a fluid phase with coexisting density η_{f} = *N*σ^{3}/*V* = 0.45 and a vacancy-rich simple cubic crystal structure with coexisting density η_{m} = 0.50 and net vacancy concentration α = (*N*_{L} - *N*)/*N*_{L} ≃ 0.064. The pressure and chemical potential at coexistence are β*p*σ^{3} = 6.16 and βμ = 18.4, respectively. The inset of Fig. 5*C* shows the equations of state and phase transitions both including and excluding the effects of vacancies in the crystal phase. The presence of vacancies in the crystal significantly lowers the melting density compared to the one reported by ref. 10, where they found that a defect-free crystal melted at η_{m} = 0.52, while the freezing packing fraction is approximately the same. We note here that we also find a melting number density of η_{m} = 0.52 if we exclude vacancies. Hence, we find that vacancies increase the range of stability of the simple cubic crystal.

### Diffusion.

As was already shown in ref. 10, the ordered phase has appreciable diffusion in the intermediate density regime. This can be understood in terms of the delocalized defects, which diffuse through the crystal and allow particles to diffuse in the opposite direction. To investigate the effect of vacancies on the diffusion coefficient in the solid, we measured the long-time self-diffusion constant of cubes in the crystal phase using EDMD simulations. Fig. 6 shows the diffusion constant as a function of density, where the net vacancy concentration at each density was chosen to correspond to the equilibrium net vacancy concentration shown in Fig. 5*C*. Near coexistence, the diffusion coefficient increases significantly, up to a maximum of *D*τ/σ^{2} = 0.05, where is the unit of time in the EDMD simulations and *m* is the mass of a single cube. For comparison, the diffusion constant in the fluid at coexistence is *D*τ/σ^{2} = 0.15, three times higher than in the coexisting solid.

At fixed density, the diffusion constant increases approximately linearly with the number of vacancies, with very little diffusion remaining at α = 0. An example of this is shown in the inset of Fig. 6, for packing fraction η = 0.56. Note that even for vanishing net vacancy concentration, diffusion is still possible via the spontaneous formation of delocalized interstitial-vacancy pairs. However, at the equilibrium net vacancy concentration (α = 0.013), the diffusion coefficient is eight times as high as in the vacancy-free crystal, indicating that vacancies play a major role in the dynamics of the particles in the solid.

## Conclusions

In this paper, we have examined the effects of vacancies on the phase diagram of hard cubes using both EDMD simulations as well as MC simulations. From the molecular dynamics simulations it is clear that vacancies play an important role in the equilibrium phase behavior of hard cubes. Free-energy calculations show conclusively a first-order phase transition between a fluid phase and a vacancy-rich simple cubic crystal phase with up to 6% vacancies. Up to the system sizes we have studied (40^{3} particles), we find long-range positional order for systems with an equilibrium concentration of vacancies (see Fig. 4). Thus, we find that the stable phase is a simple cubic crystal for all densities, albeit one with significant diffusion due to the high defect concentration.

The number of vacancies in this system is orders of magnitude larger than typically seen in colloidal systems. The stability of this vacancy-rich phase can most likely be attributed to the delocalization of defects in the crystal (Fig. 3): Clearly, one vacancy provides additional free volume for multiple nearby particles, decreasing the entropic cost of creating a defect. In most other colloidal crystals, such as hard-sphere face-centered cubic crystals, particles near a vacancy are still confined to their lattice site by their remaining neighbors. As a result, the local entropy gain from a defect is much lower. The only other similar result we are aware of is for parallel hard cubes; a somewhat artificial system that shows a very peculiar second-order freezing transition from a fluid to a simple cubic crystal (15⇓–17).

The existence of a vacancy-stabilized simple cubic phase in hard cubes leads to the question of whether vacancy-stabilized crystal structures are present in other anisotropic, entropy-driven systems. We would expect vacancies to be relevant for other (likely hard) systems with crystal structures where vacancies can delocalize. Note that delocalization also requires the absence of strong interactions that constrain particles to their lattice sites. For example, we would not expect high vacancy concentrations to occur in the simple cubic structure studied by Rechtsman et al. (18), which resulted from isotropic interactions. Recently, the phase behavior of a large number of polyhedral shapes has been studied using MC simulations. (10, 19, 20) Because vacancies are easily overlooked in the case of spontaneously formed crystals, and unlikely to form in simulations starting from a fully filled lattice, it is possible that high equilibrium vacancy concentrations occur in many of these systems. Specifically, for crystal structures where some or all of the neighboring particles can freely move into an empty lattice site, the possibility of crystal vacancies should likely be taken into account. Examples include the crystal structures predicted for hexagonal and triangular prisms in ref. 10, or the 2D (rounded) hard squares studied in refs. 21 and 22.

## Materials and Methods

### EDMD Simulations.

Please refer to *SI Text* for a full description.

### Free Energy of the Liquid.

Thermodynamic integration allows one to calculate the free energy for all densities assuming that both the equation of state and the free energy at a reference density are known. When the free energy of a reference density *F*(ρ_{0}) is known, the free energy as a function of number density *F*(ρ) can be determined using the equation of state. In particular, the free energy is given by [4]where ρ is the density and β = 1/*k*_{B}*T* with *k*_{B} Boltzmann’s constant and *T* the temperature. To measure the free energy of the fluid at a reference density, we used Widom insertion test particle method (13). The free energy of the fluid at density ρ_{0} is then given by [5]

### Solid Free Energies With and Without Vacancies.

To calculate the Helmholtz free energy as a function of the density for the solid phase, we use thermodynamic integration (13) in MC simulations of systems with *N*_{L} = 20^{3} = 8,000 or *N*_{L} = 25 × 20 × 18 = 9,000 lattice sites. We checked during our simulations that the number of lattice sites did not change spontaneously. For systems with the same density and net vacancy concentration, the differences in free energy between these two lattice sizes were within the error of our measurements. However, while equivalent free-energy calculations for a smaller system (*N*_{L} = 1,000 lattice sites) yielded qualitatively similar results, finite-size effects were noticeable when compared to the larger systems.

For the reference free energy of a crystal without vacancies, we use a variation on the method introduced by Frenkel and Ladd (23), where particles are tied to their respective lattice sites with springs, transforming the crystal into a noninteracting Einstein crystal for a sufficiently high spring constant λ. In this case, we also add an aligning potential to handle the orientational degrees of freedom of the particles (14). Using the same coupling constant λ that attaches the particles to their lattice sites, the aligning potential is given by [6]where is a unit vector along the *x*(*y*) axis, and **u**_{i,j}, with *j* = 1,2,3 are three mutually perpendicular face normals associated with particle *i*. Also, β = 1/*k*_{B}*T* is the inverse thermal energy, where *k*_{B} is Boltzmann’s constant and *T* the temperature. The parameter λ controls the strength of the external potentials; hence, for λ = 0 the system reduces to pure hard cubes, and for λ = λ_{m} with λ_{m} sufficiently large, the particles in the crystal are noninteracting.

To calculate the free energy of a system with vacancies, instead of fixing the particles to a specific lattice site, we attach the particles to their nearest lattice site (24) using [7]where **r**^{0}(**r**_{i}) is the position of the lattice site nearest to **r**_{i}. In this case, the dimensionless free energy per particle (*f* = β*F*/*N*) of the noninteracting system is , where the first term is the translational free energy of a normal Einstein crystal (13), the second term is the rotational free energy of the crystal (14), and the third term is the combinatorial entropy associated with placing *N* particles on *N*_{L} lattice sites: *f*_{comb} = -(1/*N*) log[*N*_{L} !/(*N* !(*N*_{L} - *N*) !)]. The full free energy of the crystal of hard cubes with vacancies is then given by [8]In contrast to the free energy calculations for systems without vacancies (13), the center of mass of the system is not fixed in these simulations. To equilibrate the position of the center of mass, we introduce MC moves that collectively translate every particle in the system (24). Additionally, moves that translate a single particle by exactly one lattice vector are introduced in order to improve sampling of different distributions of vacancies over the crystal. For a system with full lattice site occupancy (*N* = *N*_{L}) and thus no vacancies, we obtain good agreement between the two methods.

### Diffusion Constants.

To measure the long-time self-diffusion constant in the crystalline phase shown in Fig. 6, we performed EDMD simulations in systems of *N*_{L} = 8,000 lattice sites, for a range of densities. The vacancy concentration was chosen to correspond to the equilibrium vacancy concentration shown in Fig. 5*C*. The diffusion constant was calculated from the slope of the mean squared displacement as a function of time. An example of a plot showing the mean squared displacement is shown in *SI Text*.

## Acknowledgments

The authors thank M. Miller and D. Frenkel for many useful discussions. L.F. acknowledges support from the Engineering and Physical Sciences Research Council, United Kingdom, for funding (Programme Grant EP/I001352/1),F.S. and M.D. acknowledge support of a Vici grant from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek. M.M. acknowledges support of the SFB-TR6 program (project D3).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: L.C.Filion{at}uu.nl.

Author contributions: M.D. designed research; F.S., and L.F. performed research; F.S., L.F., and M.M. analyzed data; and F.S., and L.F. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. A.Z.P. is a guest editor invited by the Editorial Board.

See Commentary on page 17728.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1211784109/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Gottschalk S,
- Lin MC,
- Manocha D

- ↵
- Agarwal U,
- Escobedo F

- ↵
- ↵
- ↵
- Frenkel D,
- Smit B

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Damasceno PF,
- Engel M,
- Glotzer SC

- ↵
- ↵
- Wojciechowski KW,
- Frenkel D

- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics

### Related Articles

- Full of vacancies- Oct 17, 2012