Characterization of mono-and divacancy in fcc and hcp hard-sphere crystals

We determine and compare the thermodynamic properties of mono- and divacancies in the face-centered-cubic and hexagonal-close-packed hard-sphere crystals via a modified grand canonical ensemble. Widom-type particle insertion was employed to estimate the free energy of formation of mono- and divacancies, and the results are supported by an alternative approach, which quantifies the entropy gain of the neighbor particles. In hcp crystal, we found a strong anisotropy in the orientational distribution of vacancies and observe an eightfold increase in the number of divacancies in the hexagonal plane compared to the one in the out of plane at highest density of interest. This phenomenon is induced by the different arrangement and behavior of the shared nearest neighbor particles, which are located at the same distance from each vacant site in divacancy. The effect of divacancies on the free energy is to reduce that of the hcp crystal relative to the fcc by around 7 x 10(-6)kBT at melting.


I. INTRODUCTION
There exist infinitely many distinct close-packed hardsphere ͑CPHS͒ crystals that differ only in the difference in staking of the hexagonal layers.The simplest distinct stackings are the face-centered-cubic ͑fcc͒ and hexagonal-closepacked ͑hcp͒ lattices, which exhibit ABCABC... and ABAB. . .stackings, respectively.The local structural similarity of fcc and hcp crystals results in thermodynamic properties that are very similar, with the exception of some elastic constants. 13][4][5] However, there appears to be no simple theory, which can account for the numerical results. 6,7The early numerical studies in comparison of both phases have considered only defect-free systems, which excludes the fact that the inclusion of monovacancy can vary the stability of the perfect fcc crystal. 8 recent study has shown that the relative thermodynamic stability of fcc and hcp hard-sphere crystals is not affected by the inclusion of monovacancies. 9The average concentrations of divacancies in hard-sphere crystals are exponentially smaller than those of monovacancies, 10 hence the contribution of divacancies to the overall Helmholtz potential is expected to be negligible.However, in experiments, colloidal crystals often have a relatively high, nonequilibrium concentration of vacancies.Depending on their interaction, these vacancies may aggregate into specific polyvacancy structures from a lowest order, which is divacancy.In this context, it becomes important to know the formation free energies of divacancy in different structures.In particular, a strong anisotropy exists in the orientational distribution of vacancies in CPHS crystals.The approximate cell cluster theory has been used to determine the free energies and concentrations of mono-and divacancies of the CPHS crystals. 11he approached used in these theories, however, only becomes exact at close packing, which is inadequate for understanding the behavior of vacancy structures away from the close-packed limit.So far, relevant theories, which specifically deal with divacancy, ͓i.e., in accordance with arrangement and movement of the nearest neighbor ͑NN͒ particles͔ have not been found.
In the present study, we address the aforementioned issues through developing a theoretical path, which follows the grand canonical ensemble.We use molecular simulation methodologies to study mono-and divacancies in equilibrated CPHS systems; the conventional collision-based molecular dynamics ͑MD͒ method 12 is used in conjunction with the insertion-deletion bias method introduced by Pronk and Frenkel. 8In Sec.II, we elucidate the theoretical approach of the grand canonical scheme, specifically dealing with divacancy in the viewpoint of the introduction of a NN vacancy in the presence of monovacancy.A new approach to predict the free energies of divacancies is introduced.In Sec.III, we describe computational details.In Sec.IV we present and discuss results in terms of mono-and divacancies, followed by summary in Sec.V.

A. Concentration of divacancy
A statistical mechanical analysis has revealed a grand canonical potential modified to adopt a notion of the number lattice and its fluctuation. 8Without loss of generality, an extended formalism can be made for equilibrated NN vacancies with the presence of monovacancies in a system.Independency among vacancies ͑i.e., mono-and mono-, mono-and di-, and di-and divacancies͒ is assumed due to their existence in small numbers.The total free energy F 12 of a system containing monovacancies and NN vacancies is defined as , where f 0 is the free energy of a particle in a perfect crystal, −f 1 is the free energy of formation of a monovacancy at a lattice site, −f 2 is the free energy of formation of a NN vacancy, n and m are the numbers of monovacancies and NN vacancies, respectively, and M is the total number of lattices.The last term in the expression is to account for the change of the m number of the existing monovacancies into the same number of NN vacancies since the introduction of a NN vacancy triggers the change of an existing monovacancy into a NN vacancy.We define F 12 as F n − m͑f 1 −2f 2 ͒, where F n is the free energy of a system including monovacancies and independently treated throughout the derivation.
The grand canonical potential is then expressed as follows: where MЈ is M − n, ␤ is the inverse temperature 1 / k B T ͑e.g., unity in this work͒ with Boltzmann's constant k B and the temperature T , is the chemical potential, and the last term is a canonical partition function representing an approximate number of ways, by which the NN vacancies can be distributed over MЈ.The maximum of Eq. ͑1͒ appears from summing over all the states, which ⌶ M Ј visits with variation of the number of lattices ⌬MЈ.The thermodynamic relation reveals the average number of the NN vacancies ͓i.e., ignoring the terms in ϳO͑͑m / MЈ͒ 2 ͔͒ as follows: Equation ͑2͒ is valid provided monovacancies exist within the system.Thus, the product of the concentration of monovacancies ͓i.e., x MV = exp͓−␤͑ − f 1 ͔͔͒ 8 and m / MЈ becomes the average number of divacancies x DV as follows: It is important to notice that Eq. ͑3͒, which does not include −f 1 , infers that single divacancy can be permitted into the system in the course of the simulation to obtain −f 2 .The formation free energy of divacancy is estimated from twice of −f 2 .

B. Formation free energies of mono-and nearest neighbor vacancies
The formation free energy of monovacancy can be estimated by Pronk and Frenkel's equation that collects Boltzmann's factor of the potential energy change, which is sampled over the phase space of the system containing monovacancy. 8,13One can modify the formalism to imitate the simulation methodology applicable to the system con-taining divacancy ͑see Sec.III͒.Alternatively, we can estimate the formation free energies of mono-and divacancies via thermodynamic integrations using the formation entropies of those vacancies.Assuming that hard-sphere particles are confined in the Wigner-Seitz ͑WS͒ cell, then the degeneracy, in which a particle can occupy the small region within the WS cell, becomes 0 / WS , where 0 is the volume accessible by the center of the particle where the subscript "0" is denoted as the perfect crystal and WS is the volume of the WS cell.The configurational entropy of a particle is then expressed as follows: Since 0 is restricted by the movement of the particle against the NN particles, the removal of the particle induces the variation of 0 but the varying amount is in question.Note that the effective shape of 0 is assumed to be spherical.A direct quantification can be done by summing over the changes of the configurational entropies of all the surrounding particles.Conveniently, we can obtain the limit of 0 upon removal of the particle by finding the accessible volume of the center of the insertion particle.We denote the configurational entropy as s 1 , which is estimated by where 1 is the accessible volume of the insertion particle.
The difference of Eqs.͑5͒ and ͑4͒ reveals the formation entropy of monovacancy ⌬s MV , which is k B ln͑ 1 WS / 0 ͒, thus the formation free energy of monovacancy becomes Equation ͑6͒ illustrates when a particle is removed, the system reduces its energy by the free energy of that particle and the excess entropy energy.When a NN particle is removed from its site around monovacancy, which undergoes the transition to be a NN vacancy, the accessible volume 1 varies.The difference is then treated as the formation entropy of the NN vacancy ⌬s NV , which is defined as k B ln͑v 2 / v 1 ͒.Thus, the formation free energy of a NN vacancy can be expressed as follows: . ͑7͒

III. SIMULATION DETAILS
As a means to evaluate −f 2 directly, we apply Pronk and Frenkel's biased Widom particle insertion methodology 8 to a NN vacancy but with different mechanism of the simulation.The formalism is similar as follows: where ͗e −␤⌬U ͘ is the average of Boltzmann's factor of the potential energy change, which is sampled over a phase space of the system containing divacancy, ⌳ 3 is the de Bro-glie wavelength in the three dimensions, P acc ͑v WS ͒ and P acc ͑v͒ are the acceptance probabilities of the trial insertion of the fictitious particle into v WS and a biased subvolume v, respectively, and P rem ͑v͒ is the acceptance probability of a particle resides in v.The first expression in Eq. ͑8͒ can be simply interpreted as the accessible volume in the WS cell as shown in the second expression.It is inefficient to choose the entire region of the WS cell for the insertion trial, thus the subvolume is introduced and this scheme requires another quantity, which compensates the proportional difference due from the reduction of the insertion region.
To obtain P acc ͑v͒, we introduce a divacancy in the system and perform the insertion trial into a NN vacancy.For P rem ͑v͒, a monovacancy is introduced and the probability of a NN particle, which happens to be in v, is collected.Note that the probability represents the fractional time since we employ MD simulation.We call this approach the insertiondeletion simulation ͑IDS͒ hereafter.The spherical volume of the 80% radius of a sphere inscribed into the WS cell is adopted in this work.To fix the orientation of divacancy, the single occupancy cell ͑SOC͒ method is used to confine particles with the diameter of the inscribed sphere.Our preliminary study revealed that the average mean square displacement of the particles does not vary when the confining diameter is set to be above 90% of the NN distance.Variation of density is of our interest to observe the trend of the behavior of the NN particles and its effects on the properties of vacancies.The range of the density is from * = 1.2 to melting, 1.0376, where * is M 3 / V with the total volume of the system V, the diameter of a hard-sphere ͑e.g., unity in this work͒, and M is the number of lattices in the simulations.We use 576 ͑6 ϫ 8 ϫ 12 unit cells͒ to 1728 ͑12ϫ 12ϫ 12 unit cells͒ lattices, which are constructed in rhombic boxes to avoid size effects and to accommodate the correct behavior of the NN particles.The chemical potential and f 0 ͑see Fig. 1͒ is determined by thermodynamic integration using Speedy's equation of state 14 from the melting point provided by Frenkel and Ladd. 2

A. Monovacancy
The estimates from Eq. ͑6͒ agree well with f 1 's independently obtained by IDS approach and the values of ⌬s MV are reported in Table I.As density increases, the rate of the entropy gain by the neighbor particles reduces but the removed particle possessed higher free energy in the system so that −f 1 decreases.Figure 1 shows the density dependence of −f 1 , which is obtained from IDS approach.In comparison of both fcc and hcp crystals, there exist little discrepancy in −f 1 , separately the calculated values of ⌬s MV are approximately equal within confidence limits, thus structural effects of the NN particles on monovacancies are indistinguishable.We illustrate this phenomena under assumption that the NN particles out of entire neighbors gain most of entropy ͑e.g., about 80%͒ 15 by removal of a particle.The average distribution function P͑r * ͒ of the NN particles ͑e.g., without SOC͒ is constructed by weighting the collection of the distance of a NN particle ͑i.e., the fractional time spent in the predeter-mined shell volume by the NN particle͒ from its original lattice site with spherical shell volumes.As shown in Fig. 2, close overlap of P͑r * ͒'s for the fcc and hcp crystals shares the same overview done by comparing ⌬s MV 's; little effects of anisotropic arrangements of the NN particles on monovacancies are found.
In comparison of P͑r * ͒'s at different densities in vacancy-free crystals with those in crystal with monovacancy, we found abrupt change of P͑r * ͒, which indicates the different behavior of the NN particles depending on density.In particular, the significant variation is observed over the range of the density from 1.1 to 1.2.By projecting the relative positions of the NN particles from their original lattices into each directional vector connecting the lattice sites with the vacant site, the existence of uneven arrangements of the NN particles was found.At * = 1.17 for the fcc crystal and 1.16 for the hcp, 6 ϫ 6 competing ͑e.g., number indicates the number of the NN particles͒ movement of the NN particles is initiated followed by the 2 ϫ 4 ϫ 6 competition beginning at * = 1.15 and ending at 1.1.Reduction of density slowly provides limiting spaces for some of the NN particles to be moved more freely than the others around monovacancy.Below * = 1.1, those unique competitions diminish.From independent simulations, we observed that the particles participated in those specific movements are random and the uneven arrangements of particles are hardly varied over a simulation.Despite different arrangements and movements of each NN particle ͓i.e., different P͑r * ͒'s for each NN par-ticle͔, the averaged P͑r * ͒'s at different densities correspond to the order of ⌬s MV from * = 1.1 to 1.2.Pair correlation function of the surrounding particles ͑work not shown͒ illustrates the relaxation of the particles beyond the NN particles FIG. 1. Plot of the free energies vs *.Open and filled markers represent the fcc and hcp systems, respectively.The different markers represent the different types of free energies; circle, triangle, square, and diamond represent −f 0 , −f 1 , and −f 2 of the NN vacancies in divacancy in the fcc crystal and divacancy oriented in adjacent planes ͑i.e., AB plane͒ of the hcp, and −f 2 of the NN vacancy in divacancy oriented in the same plane ͑i.e., AA or BB plane͒ of the hcp, respectively.Lines are a guide to the eye.The 67% confidence limits are smaller than the size of the markers.so that the range of P͑r * ͒ becomes narrower at the density close to melting in comparison of the one at * = 1.1.We have not observed the same unique change of P͑r * ͒ for the case of divacancy in both crystals.This is due to excess space brought by introducing the NN vacancy around monovacancy, which reduces the competitions of the NN particles, thus the density around 1.2 is not high enough to trigger the similar competitions as in the case of monovacancy.

B. Nearest neighbor vacancy in divacancy
With variation of density, we observed that two distinct values of −f 2 's arise within the hcp crystal, as shown in Fig. 1.The formation free energies of a NN vacancy from both methods from Eq. ͑7͒ and IDS approach agree well each other, and ⌬s NV is reported in Table I.In particular, concurrence of −f 2 's in fcc and hcp crystals postulates the existence of the structural similarity of some of the NN particles.Out of 18 NN particles, the four particles shared by the NN sites in divacancy ͑i.e., the particles, which are at the same distance from each site in divacancy͒ play an important role to induce the similar property.If divacancy is oriented through adjacent planes ͑i.e., AB͒, the lattice sites of the shared particles are positioned on the plane perpendicular to the divacancy orientation and depict a rectangular shape.We found that this shape is exactly the same in the fcc crystal at a given density.On the other hand, for the orientation of divacancy in a same plane ͑i.e., AA or BB͒, the lattice sites of the shared particles make a kite shape, which can be projected on a hexagon; the three particles locate at the head side against one, which locates at the tail side.
We resort to P͑r * ͒ to illustrate the behavior of the binding strength.For divacancy in the fcc crystal at * = 1.2 as shown in Fig. 3͑a͒, the shifted centers of P͑r * ͒'s to the right indicates that the seven nonshared NN particles spend more time at the region ͑i.e., r * Ϸ 0.2͒ toward a NN vacancy in divacancy with the four shared ones toward the center of the divacancy.The behavior of P͑r * ͒'s illustrates the tight wrapping of the each 11 NN particles on the each NN vacancy, which affects to the reduction of the interactions of the NN vacancies in divacancy at high density.At lower density * = 1.05, the centers of P͑r * ͒'s become around the original lattice sites of the NN particles due to the relaxation of the surrounding particles beyond the NN ones.The tails of P͑r * ͒'s become short ranged in comparison of the ones at high density.This phenomenon consequently induces the more interactions of the NN vacancies in divacancy.Similar behavior of the NN particles is observed for the divacancy in the AB plane of the hcp crystal, as shown in Fig. 3͑b͒.While the NN particles do not gather around divacancy at high density, the NN particles move wider range and the tails of them match with the ones in the case of the fcc.At low density, P͑r * ͒'s from both fcc and hcp cases overlap with each other.Thus at both densities, we can expect the similar binding strength of divacancy.The associated free energy f a ͓e.g., estimated by 2͑f 2 − f 1 ͔͒ indeed shows weak binding of  divacancy at high density and becomes stronger at low density, and the values are reported in Table I.When divacancy is oriented in a plane of the hcp crystal, f a tends to increase as the density increases.P͑r * ͒'s appeared to be similar to those shown in Fig. 3͑b͒ except one for the particular shared particle, as shown in Fig. 3͑c͒.The particle was initially located at the tail side of a kite shape.While the other NN particles wrap the divacancy with a same degree as in the case for the divacancy in the AB plane, the particular particle penetrates deeper into the center of divacancy and induce the increase of f 2 .This phenomenon becomes severe at high density, thus higher value of f a is obtained.The existence of mono-or divacancies reduces the free energy of the crystal as density increases.However, the net differences of the chemical potential and the formation free energies of the vacancies ͓i.e., − f 1 or − f 2 , which was defined as the driving force to insert particle͑s͒ to remove vacant site͑s͒ 13 ͔ offset the increasing stabilities of the vacancies as density increases, thus the vacancies result in small numbers, which are reported in Table I.We found manifold increase of the number of divacancies in the AA ͑or BB͒ plane of the hcp crystal by ϳ160% and ϳ790% at * = 1.05 and 1.2, respectively, compared to the ones in the AB plane.Our estimates agree well with available values in literature except Langeberg et al.'s, 11 the extension from the close-packed limit leads large deviations especially within the hcp crystals, indicating obsoleteness of two-particle cell cluster method at low density.At melting, the existence of divacancy shifts the chemical potential by around 1.3 ϫ 10 −5 k B T and 2.0ϫ 10 −5 k B T for divacancy in the AB plane ͑i.e., similar in the fcc crystal͒ and the AA ͑or BB͒ in the hcp at melting, respectively.In comparison to the case of monovacancy ͑i.e., around 2.3ϫ 10 −3 k B T in both crystals͒, 9 the existence of divacancies may not affect to alter melting point and relative stability of the fcc and hcp crystals.
While mobility of vacancy is categorized as rare event, it is worth noting that the behavior of the NN particles is closely involved in diffusion and dissociation events.By examining the tails of P͑r * ͒'s, there exist small yet apparent differences in monovacancy-mobility favoring in especially at * = 1.1, as shown from Fig. 2. As to divacancy, dissociation, and diffusion compete with ratio of 3.5 to 1.For divacancy in the fcc crystal, the dissociation of divacancy requires the free energy of around 0.45k B T ͓i.e., f 2 − f 2 ͑r second ͒, where r second is the distance from a vacant site to the second neighbor vacant site͔.Continuous dissociations of the separated neighbor vacancies are unfavorable due to raising the free energy of the system, but the thermal energy 16 is still compatible for particles to participate in separating the neighbor vacancies.For divacancy in the AB plane of the hcp crystal, there exist two resulting orientations of divacancy by different hopping mechanisms of the four shared particles; from the AB plane to the AB or to the AA ͑or BB͒.Pairs of diagonal particles in the shared ones ͑i.e., the rectangular shape͒ show same probability of hopping of the particles into divacancy, thus resulting orientation has same chance to be in the AB plane or the AA ͑or BB͒.We found that the diffusion of divacancy in the path of the AB to the AA ͑or BB͒ is favored than the one in the reverse path, which requires more free energies to overcome the difference in the association free energies ͑i.e., ϳ2.0k B T and ϳ0.5k B T at * = 1.2 and 1.05, respectively͒.For divacancy in the AA ͑or BB͒ plane, dissociation of divacancy occurs less than diffusion, which particularly happens within that plane due to the penetrating particle located on the same plane.The other three shared NN particles participate less in diffusion due to colliding and hindering competitions among themselves in accordance with the nonshared NN ones.

V. CONCLUSION
Modification of Pronk and Frenkel's grand canonical ensemble enables one to obtain the free energy of formation of a NN vacancy in divacancy.Our analysis is focused on comparing the thermodynamic properties of mono-and divacancies in accordance with analyzing the behavior of the NN particles.We developed an alternative path, which uses the configurational entropy to calculate the formation free energies of mono-and divacancies and found a good agreement in results with IDS approach.There exist little overall structural effects of the 12 NN particles on monovacancy in the fcc and hcp crystals over the density of interest.At high density, the gathering of the NN particles induces lower −f 1 .As density decreases, the NN particles uniquely relax from * = ϳ 1.17 to 1.1, and that unique movement of the particles diminishes at close to melting due to relaxation of the sur- rounding particles.For divacancy, the effects of the arrangements of the four shared NN particles become significant on thermodynamic properties of divacancy.In particular, we found that the behavior of the shared ones induces the considerable discrepancy within the hcp crystal.The shared NN particles depict the shapes of rectangle ͑i.e., same for the shared ones in the fcc lattice͒ and kite on hexagon, which correspond to the arrangements of divacancies in the AB plane and the AA ͑or BB͒, respectively.The formation free energy of divacancy in the later case found to be lower due to wide movement of a shared particle, which results in the increase of the number of divacancies.In hcp crystal, we estimated an eightfold increase in the number of divacancies in the AA ͑or BB͒ plane compared to the one in the AB at * = 1.2.This result negates Langeberg et al. prediction, which shows higher concentration of divacancies oriented in the AB plane, extrapolated from the close-packed limit.In comparison of dissociation and diffusion events, analysis of thermodynamic stability shows preference in the diffusion in both crystals.Within the hcp lattice, the diffusion in the same plane is favored.
The analysis presented in this work could be applied to any model systems provided the chemical potential and the equation of state are known.The developed ensemble could also serve as a seed formalism to study specific polyvacancies ͑e.g., especially in a form of chain, which contains less complicated interactions among vacancies͒ via modifying the exponent in Eq. ͑3͒.Implementation of the formalism could be rather straightforward but must include the statistical mechanical information of lower order vacancies.Recently, Torquato and Stillinger 17 investigated tunneled ͑i.e., chain vacancy͒ crystals on packing threshold and showed the existence of a new hard-sphere crystal structure at low density with close-packed limit.The chain type of polyvacancy formed from self-diffusion is theoretically possible, yet the cluster type is favorable due to lowering the free energy of the crystal 10,16 with respect to the equilibrated number of vacancies.The mechanically stable tunneled crystal contains extremely large quantity of vacancies beyond equilibrium.To compare with thermodynamic stability, one must consider the thermodynamic limit of excess number of vacancies, which is still in an open question.

4 2 ϫ 10
FIG. 2. Average distribution functionsP͑r*͒ of the NN particles centered on their lattice sites around monovacancy, where r* is r / .Open and filled markers represent the fcc and hcp systems, respectively.The different types of markers represent different densities; circle, square, triangle, and inverse triangles are assigned as * = 1.2, 1.1, 1.05, and 1.0376, respectively.An upper-right plot expands P͑r*͒ at * = 1.2.Lines are a guide to the eye.The 67% confidence limits are smaller than the size of the markers.

FIG. 3 .
FIG. 3. Average distribution functionsP͑r*͒ of the NN particles centered on their lattice sites around a NN vacancy, where r* is r / .Diamond and circle markers represent * = 1.2 and 1.05, respectively.For the plot ͑a͒, open and filled symbols represent the seven NN particles and the four shared NN ones, respectively.For the plot ͑b͒, the same notations are applied and divacancy is oriented in the adjacent planes ͑i.e., AB plane͒ of the hcp crystal.For the plot ͑c͒ divacancy is oriented in the same plane ͑i.e., AA or BB plane͒ of the hcp crystal.The same notations are applied except the filled markers representing the three shared NN particles ͑i.e., on head side of a kite shape͒.For the one shared particle ͑i.e., on tail side of a kite shape͒, open triangle and inverse triangle were used for * = 1.2 and 1.05, respectively.Lines are guides to the eye.The 67% confidence limits are smaller than the size of the markers.

TABLE I .
Concentrations and association free energies of divacancies, and formation entropies of mono-and NN vacancies in the fcc and hcp CPHS systems.͑Reported values are based on around 4 ϫ 10 4 collisions per particle.Subscripted digit represents the statistical error at the last digit.͒