# Molecular Dynamics Simulation of Binary Hard-Sphere Crystal/Melt Interfaces

## I Introduction

In order to fully understand such important phenomena as near-equilibrium crystal growth, homogeneous nucleation and interfacial solute segregation, a detailed microscopic description of the crystal/melt interface is necessary. Since direct experimental data is scarce, due to the extreme difficulty of constructing and interpreting experiments on such systems, computer simulation has become an important tool, not only in its usual role of aiding in the development and evaluation of interface theories, but also in determining the basic microscopic phenomenology of crystal/melt interfaces. Previous simulation studies of such interfaces have focused entirely on single-component systems (for a review of recent simulation studies, see Ref. ); however, since most technologically important materials are not pure substances but mixtures (such as alloys), it is important that such studies be extended to multi-component interfacial systems. In addition to the usual issues of interfacial width, interfacial free energy, and transport within the interface, multi-component systems allow for the study of interfacial segregation.

As in single-component systems, the crystal/melt interface of a multi-component system is characterized by measuring the change in the various structural, thermal and dynamical properties of interest as one traverses the interface from one bulk phase into the other. For planar interfaces (the type studied here), the $`z`$ axis is usually taken as the direction perpendicular to the interfacial plane and quantities are averaged over $`x`$ and $`y`$ and presented as functions of $`z`$. Examples include the density profiles of the various components (labeled by $`i`$), $`\rho _i(z)=\rho _i(𝐫)_{xy}`$ and the diffusion constant profiles, $`D_i(z)`$. The thermodynamic quantity of most interest is the solid-liquid interfacial free energy, $`\gamma _{\mathrm{sl}}`$, which is defined as the work required to form one unit area of surface – this quantity is extremely difficult to calculate via simulation – even in the single-component case. The only reliable calculation is the calculation of $`\gamma _{\mathrm{sl}}`$ for a single-component Lennard-Jones system by Broughton and Gilmer.   

Recently, we have reported initial simulation results on the crystal/melt interface of a two-component hard-sphere mixture. Specifically, the interface between the (100) face of a disordered face-centered cubic (FCC) crystal and the coexisting fluid was studied. In this work we extend this calculation to the (111) face and revisit the (100) interface calculations using more detailed analysis. These simulations represent the first such simulations on the crystal/melt interface of a multi-component system.

The reasons for choosing the hard-sphere system for this initial study are two-fold: First, it is now well established that the structure and freezing behavior of dense, simple fluids, is determined, for the most part, by packing considerations determined by the repulsive part of the interaction potential. The effect of the attractive forces can generally be accounted for by treating it as a perturbation to the repulsive part of the potential, which is often approximated by a hard-sphere with some effective diameter. Thus, the hard-sphere system is a useful reference from which to begin studies of more realistic systems. Studying the hard-sphere system also allows one to directly probe the role of packing (which is purely entropic) in determining the interfacial phenomenology. Second, the relative simplicity of the hard-sphere system lends itself well to theoretical study – the vast majority of density-functional calculations on crystal/melt interfaces involve hard-sphere systems. The disordered-FCC/fluid interfaces were chosen for initial study, because, in order to begin an interface simulation, the phase coexistence conditions must be very accurately known, and the disordered-FCC/fluid region of the binary hard-sphere phase diagram (which occurs in the region of the phase diagram where the difference between the diameters of the different components is not too great) has been well characterized.

In the next section we describe the binary hard-sphere system and its phase behavior. The results of the disordered FCC/fluid interfacial simulations reported for the (100) and (111) interfaces are discussed in Section 3. In Section 4, we conclude.

## II The Binary Hard-Sphere System

We consider a two-component system consisting of hard spheres of differing sizes. The interaction potential for such a system can be written

$$\varphi _{ij}(r)=\{\begin{array}{cc}\mathrm{}\hfill & r<\sigma _{ij}\hfill \\ 0\hfill & r\sigma _{ij}\hfill \end{array},$$
(1)   
where $`r`$ is the distance between the centers of two spheres and $`i`$ and $`j\{1,2\}`$ index the two types of spheres, which are distinguished by their different diameters, $`\sigma _1`$ and $`\sigma _2`$, while their masses are assumed to be identical. We also assume that the spheres are additive; i.e., $`\sigma _{11}=\sigma _1`$, $`\sigma _{22}=\sigma _2`$, and $`\sigma _{12}=\sigma _{21}=(\sigma _1+\sigma _2)/2`$. The following definitions and conventions are adopted for this study. First, it is assumed, without loss of generality, that $`\sigma _2\sigma _1`$, and the hard-sphere diameter ratio $`\alpha `$ is defined as

$$\alpha =\sigma _1/\sigma _2,(0\alpha 1).$$
(2)   
Second, if there are $`N_1`$ hard spheres with diameter $`\sigma _1`$ and $`N_2`$ with diameter $`\sigma _2`$ in the volume $`V`$, then

$$\rho =\frac{N_1+N_2}{V}=\rho _1+\rho _2$$
(3)   
is the total number density, and $`\rho _i`$’s represent the respective number densities for individual species. Third, the concentrations (mole fractions) of individual components are given by

$$x_i=\rho _i/\rho ,i=1,2.$$
(4)   
Since $`x_1+x_2=1`$, a single variable $`x`$ is usually used, such that $`x_2=x`$ and $`x_1=1x`$.

Also, the total packing fraction for the mixture $`\eta `$ is

$$\eta =\eta _1+\eta _2,$$
(5)   
where $`\eta _i=\frac{\pi }{6}\sigma _i^3\rho _i`$, $`i=1,2`$. The diameter ratio $`\alpha `$ together with any pair of independent parameters from those defined above can be used to completely specify the fluid state of a binary hard-sphere system. The unit of length for the binary system is taken to be equal to the diameter of a larger sphere $`\sigma _2`$.

Depending on the value of the diameter ratio $`\alpha `$, different crystal structures may have the lowest free energy . In particular, for $`\alpha `$ above about 0.85, the substitutionally disordered FCC crystal, in which spheres of different diameters are distributed randomly over the sites of an FCC lattice, is the stable structure at freezing. For $`\alpha `$ below this value, a variety of structures may exist including ordered solid such as NaCl and CsCl structures. For an interface simulation study, it is necessary that the phase coexistence conditions be accurately known. Therefore, we have chosen for this study the disordered-FCC/melt system, since the coexistence conditions for the equilibrium disordered crystal/melt interface can be found in the study of Kranendonk and Frenkel, who have calculated the crystal/melt phase diagrams for several values of the diameter ratio in the range $`0.85<\alpha <1`$. For the present study, the diameter ratio $`\alpha =0.9`$ is selected.

For a binary system with a given diameter ratio $`\alpha `$, the coexistence state is specified by two number densities $`\rho ^\mathrm{c}`$ and $`\rho ^\mathrm{f}`$, as well as two concentrations $`x^\mathrm{c}`$ and $`x^\mathrm{f}`$ (or, alternatively, by four number densities of individual species, $`\rho _i^\mathrm{c}`$ and $`\rho _i^\mathrm{f}`$, where $`i=1,2`$). The pressure-concentration phase diagram determined in the simulation by Kranendonk and Frenkel for $`\alpha =0.9`$ is shown in Fig. 1. At a given pressure, the fluid and crystal phases have different concentrations at coexistence, represented in the plot by triangles and circles, respectively. To maximize the deviation from single-component behavior, a point on the $`\alpha =0.9`$ phase diagram was chosen for the interface simulation, where the concentration difference between crystal and melt phases is the largest. This point occurs at a pressure of

$$P=14.7k_\mathrm{B}T/\sigma _2^3$$
(6)   
and concentrations of

$$x^\mathrm{c}=0.71\text{and}x^\mathrm{f}=0.54$$
(7)   
for the crystal and fluid phase, respectively, and is shown in Fig. 1 by the solid triangle and circle. We have found the total packing fractions (number densities) in each two phases at

$`\eta ^\mathrm{c}`$ $`=`$ $`0.552(\rho ^\mathrm{c}\sigma _2^3=1.144),`$ (8)   
$`\eta ^\mathrm{f}`$ $`=`$ $`0.502(\rho ^\mathrm{f}\sigma _2^3=1.096)`$ (9)   
and have used these values for the initialization of the binary interface systems. We have run simulations for several trial systems and have found that the systems are stable and the bulk crystal remains on average stress-free. Therefore, the above parameters have been used for all the simulation runs from which average properties of the interface have been computed.

## III Interface Construction and Simulation Results

To create the initial bulk systems that are placed together to form the interface, the two sphere types are distributed randomly according to the crystal and fluid coexistence concentrations. The concentration in each crystal layer is maintained fixed by randomly distributing the spheres of different types on a layer by layer basis, thereby removing layer-to-layer concentration fluctuations due to finite system size. This constraint is not expected to affect the results in any significant way besides removing the fluctuations that cannot be averaged over during the simulation run due to practically no diffusion in the bulk crystal. The random distribution of particle types is justified by the conclusion of Kranendonk and Frenkel that above about $`x=0.6`$ there is little or no local substitutional ordering.

Since the hard-sphere system evolves on a collision-by-collision basis, the natural unit of time for hard spheres is the mean collision time $`\tau _\mathrm{c}`$, i.e. the average time between collisions suffered by a given particle. On the other hand, the duration of a simulation run is most conveniently measured in terms of the total number of collisions. In the present study, in order to have better correspondence between the simulation time and the system evolution time, we measure simulation time in units of the number of collisions per particle (cpp), defined as twice the ratio of the total number of collisions to the number of particles in the system, so that $`1\mathrm{c}pp\tau _\mathrm{c}`$.

We have prepared 10 systems for each of (100) and (111) crystal orientations and have run the simulations for 20 000 cpp with the interfacial diagnostics being recorded every 200 cpp. The (100) systems contain 11 616 spheres and have dimensions $`L_x=L_y=16.70\sigma _2`$, and $`L_z=37.09\sigma _2`$, while the size of the (111) systems is 11 340 spheres and $`L_x=16.09\sigma _2`$, $`L_y=16.72\sigma _2`$, and $`L_z=37.60\sigma _2`$. (We have also done simulations on these systems with smaller numbers of spheres – $``$ 3000 and 6000 – and different cross-sectional shapes. The results of these simulations are within statistical error, quantitatively identical to those presented here. The larger samples have been chosen as they give better statistics, have much shorter equilibration times than the smaller samples, and exhibit smaller interfacial fluctuations.)

The concentration fluctuations in the fluid and the interfacial regions have been found to be much larger and more persistent than the density fluctuations. Also, due to a finite system size, even though the total momentum with respect to the simulation cell is set to zero, drift of the interface positions is observed. In order to avoid broadening of the interfacial profiles caused by the drift, we have selected for the final averages 12 segments of 2000 cpp in duration from each crystal orientation, such that the drift does not exceed half the distance between crystal layers (for more details on the methods of interfacial construction and equilibration used here, as well as more detailed definitions of measured (computed) profiles, see our recent work on single-component systems ).

The fine-scale profiles for the two components $`\rho _1(z)`$, $`\rho _2(z)`$ and for the total density $`\rho (z)=\rho _1(z)+\rho _2(z)`$ are shown in Figs. 2 and 3 for the (100) and (111) crystal orientations, respectively. For the total density profiles the 10-90 widths of the height of the density peaks equal to 5.3 and 5.6$`\sigma _2`$ for the (100) and (111) interfaces, respectively. The density oscillations in $`\rho _2(z)`$ and $`\rho (z)`$ dampen monotonically, while $`\rho _1(z)`$ exhibits a peculiar non-monotonic peak-height envelope, a phenomenon that has not been seen in any of the single-component system studies.

The non-monotonic behavior of the fine-scale density, $`\rho _1`$, can be explained by examining the coarse-grained (filtered) density profiles $`\overline{\rho }_1(z)`$, $`\overline{\rho }_2(z)`$, and $`\overline{\rho }(z)`$, computed using Finite Impulse Response (FIR) filter. The use of such filters for coarse-graining the density profiles is necessary when the peak-to-peak spacing of a profile is not constant through the interface – in such a case the use of uniform bins to perform coarse graining can lead to misleading results. (For a detailed description of the use of such filters in analyzing density profiles, see our recent article on the single-component hard-sphere interface .)

The filtered profiles are shown in Figs. 4 and 5 for the (100) and (111) interfaces, respectively. In both cases the individual species densities change have 10-90 widths of about 2$`\sigma _2`$, whereas the corresponding width for the total density is about 4$`\sigma _2`$. This seems strange in that the total density is the sum of the individual densities, but since the difference between the total densities on either side of the interface is an order of magnitude smaller than that for the individual densities, very small changes in $`\rho _1`$ or $`\rho _2`$ can contribute significantly to the 10-90 width of the total density, while remaining unimportant in determining the width of the individual densities. (Note that, the width of the total density is somewhat larger than the 3.2-3.3$`\sigma `$ found in the single-component case).

The rapid change in concentration over about 2$`\sigma _2`$, combined with the fact that the density oscillations on the fluid side of the interface in the fine-scale profiles persist 2-3$`\sigma _2`$ after the concentrations have relaxed, leads to the above mentioned non-monotonicity. Since the number density of the smaller spheres in the fluid \[$`\rho _1^\mathrm{f}=(1x^\mathrm{f})\rho ^\mathrm{f}=0.504\sigma _2^3`$\] is larger than the corresponding density in the crystal region \[$`\rho _1^\mathrm{c}=(1x^\mathrm{c})\rho ^\mathrm{c}=0.332\sigma _2^3`$\], the ordering of the fluid in the presence of the interface occurs in a region with higher average density than that in the bulk crystal, resulting in the higher profile peaks in the interfacial region.

Analysis of the coarse-grained densities also leads to the conclusion that there is no statistically significant equilibrium interfacial segregation (Gibbs adsorption) in either of our simulated interfacial orientations. Such segregation is quantified by $`\mathrm{\Gamma }_1`$/A, where $`\mathrm{\Gamma }_1`$ is the excess number of type 1 (small) particles (here the smaller particle is taken to be the solute), defined using a Gibbs dividing surface for which the excess number of type 2 (large) is zero, and $`A`$ is the area of the interface. The significance of this result should be taken, however, in light of the issues of chemical equilibrium discussed below.

Information about changes in the inter-layer spacing is obtained by measuring the layer separation, defined as

$$\mathrm{\Delta }z_i=\overline{z}_{i+1}\overline{z}_i,$$
(10)   
where $`\overline{z}_i`$ is the center of mass of layer $`i`$ determined from the fine-scale density profile between the adjoining density minima. We calculate $`\mathrm{\Delta }z_i`$ for the total density profile as well as for the density profiles of individual components. The results are shown in Fig. 6 with diamonds, squares and circles representing layer separation in $`\rho _1(z)`$, $`\rho _2(z)`$ and $`\rho (z)`$, respectively. As in the single-component case, the layer separation shows large layer expansion for the (100) interface and very little expansion for the (111) interface. In addition, we see significantly different behavior of the individual density profiles. For the (100) orientation the increase in the layer separation of the $`\rho _1(z)`$ profile is delayed by about two crystal layers, compared with that of $`\rho _2(z)`$. For the (111) orientation the $`\rho _1(z)`$ profile exhibits almost no change in the interlayer spacing. Evidently, at the onset of the disorder in the interfacial region, the spheres of type 2, having larger diameter, are repelled farther away from the ordered crystal layers.

The self-diffusion constant profiles are computed separately for the two particle types according to the formula

$$D_i(z)=\frac{1}{6N_i(z)}\frac{d}{dt}\underset{j=1}{\overset{N_i(z)}{}}[𝐫_j(t)𝐫_j(t_0)]^2,i=1,2,$$
(11)   
where $`N_i(z)`$ is the number of spheres of type $`i`$ between $`z\delta z/2`$ and $`z+\delta z/2`$ at time $`t=t_0`$, where $`\delta z`$ is the crystal layer spacing, and the brackets represent the average over time origins $`t_0`$. The sphere displacement was monitored on a uniform coarse scale over time $`t_{\mathrm{max}}t_0=5.5(m\sigma _2^2/k_\mathrm{B}T)^{1/2}`$. The averaging was done over 50 time origins for each of the 12 selected intervals. The average diffusion constant profiles for the (100) and (111) interfaces are shown in Fig. 7. As could have been anticipated, the self-diffusion coefficient in the bulk fluid is larger for the smaller spheres at $`0.019(k_\mathrm{B}T\sigma _2^2/m)^{1/2}`$, compared to that for the larger spheres at $`0.016(k_\mathrm{B}T\sigma _2^2/m)^{1/2}`$. The diffusion profiles change monotonicly across the interface with the 10-90 widths of about 3.3$`\sigma _2`$ for both (100) and (111) interfaces, which is intermediate between the widths for the interfacial profiles and that for the total density. Note that, since the midpoint of the diffusion profile is shifted by about 1-2$`\sigma _2`$ to the liquid side of the midpoints of the density profiles (that is, the bulk of the increase in the diffusion constant occurs after the densities have relaxed to nearly their liquid values), the actual width of the interfacial region is larger than either transport or structural properties would indicate alone.

At this point it is useful to comment on the question of chemical equilibrium in our systems. When we create the interface by placing crystal and fluid blocks next to each other and then allowing the system to evolve, we expect that after some time the system will stabilize itself and the interfaces will be in equilibrium. In our simulations we see that after the initial equilibration, both temperature (see Ref. ) and transverse pressure profiles exhibit no significant deviation in the interfacial region from their average values. This is an indication that the interface is in thermal and dynamic equilibrium with the surrounding bulk phases. On the other hand, the concentration equilibrium of individual species at the interface cannot be assumed with the same certainty. One admitted approximation that we have made in the construction of our interface is the use of a randomly substituted solid mixture. In their studies on hard-sphere binary mixtures, Kranendonk and Frenkel report no significant local substitutional order for solid solutions above 60% large spheres. Since our solid has a large-sphere mole fraction of 0.71, our assumption of random substitution in the bulk solid should be valid. However, the work of Kranendonk and Frenkel applies strictly to the bulk system and it is possible that there is some equilibrium local substitutional order that should be present on the solid side of our interface that we cannot see due to the very long relaxation times for concentration fluctuations. The closer we are to the crystal side of the interface, the less certain can we be that a particular configuration of the small and large spheres correctly represents the equilibrium concentration profile. The reason is that, unlike temperature and pressure fluctuations which are transported through the system via collisions, the concentration fluctuations are introduced when particles themselves drift through the system. This obviously requires moderate values of the diffusion coefficient which cannot be achieved in the crystal. (In order to achieve true concentration equilibrium in a binary interfacial system one would need to simulate composition fluctuations consistent with the chemical potential balance in both crystal and fluid phases and across the interface. This can be done, for example, by introducing Monte-Carlo moves into the equilibration process that allow small spheres to become large ones and vice versa, with probabilities that produce correct chemical potential profiles. We have not, however, found a consistent way of doing this in practice that both preserves the stability of the interface and gives good statistics.) This being said, the effect of the uncertainty in the degree of chemical equilibrium achieved should have little effect on most of the phenomena mentioned above, such as non-monotonicity of the $`\rho _1`$ profile and the anomalous behavior of the lattice spacing in both interfaces, since these effects are primarily due to behavior on the liquid side of the interface where the diffusion constant is large enough to ensure relaxation in concentration. The largest effect will be on the width of the concentration profiles on the solid side of the interface - the value given above should be taken as an upper bound - and the precise degree of equilibrium solute segregation (adsorption) on the solid side. The high stability of our interfaces leads us to speculate that these effects will be small, but, because of the problems outlined above, they cannot be discounted.

## IV Summary and Conclusions

We have presented detailed molecular-dynamics simulations for the (100) and (111) (disordered) FCC crystal/melt interfaces for a binary system of hard spheres. The ratio of the small hard-sphere diameters and that of the large spheres was chosen to be 0.9. This study extends our earlier preliminary work on the (100) interface for this system. The principle results of this study are summarized as follows:

1. The fine-scale density profiles for the smaller particle, $`\rho _1(z)`$, in contrast to the single-component case, exhibit a pronounced non-monotonic envelope. This behavior is not seen in either the total or large-particle density profiles. Analysis of the coarse-grained density profiles shows that this phenomenon is not due to any appreciable adsorption of the smaller particle at the interface, but is entirely due to the fact that increase in the small particle concentration occurs over a shorter length scale than the decay of the density ocsillations in the liquid.
2. As in the single-component case we see an increase in the spacing between the (100) density-profile peaks as the interface is traversed from crystal to fluid. A much less pronounced effect is seen in the (111) interface where the large-particle density-peak spacing stays mostly constant except for a maximum well on the liquid side; in contrast, the spacing for the smaller particles is constant through the interface.
3. The widths of the coarse-grained concentration profiles (calculated with a Finite Impulse Response (FIR) filter) are considerably smaller that those for the total densities (about 2$`\sigma _2`$ versus 4$`\sigma _2`$) and no significant equilibrium interfacial segregation is seen. As in our earlier single-component study, the use of FIR filters to determine coarse-grained density profiles is necessary to avoid artifacts of the binning process when the peak separation is not constant.
4. The diffusion profiles and coarse-grained density are shifted by 1 to 2$`\sigma _2`$ relative to one another. Significant diffusion begins only after the density has relaxed to nearly the bulk liquid value. Therefore, the interfacial region is wider than either of these profiles would indicate separately, and we can identify two separate sub-regions as the interface is traversed from solid to fluid, in which density relaxation or changes in transport properties are dominant, respectively.