Epitaxial nucleation of a crystal on a crystalline surface

We consider the nucleation of a crystal phase, on a crystalline surface of a different substance. Sixty years ago, Turnbull and Vonnegut predicted that a crystalline surface is best at inducing nucleation of another crystal when there is a perfect epitaxial match between the two bulk lattices. We use computer simulation to show that this is not quite right. In fact, the crystal lattice of a finite nucleus is strained from that in the bulk, and nucleation is fastest when the surface matches this strained lattice. We also find that the approach of Hillier and Ward predicts when nucleation is epitaxial.

Simulation setup. -We study two surfaces: an hcp (hexagonal close packed) lattice with an exposed (0001) plane, and an fcc (face centered cubic) lattice with an exposed (100) plane. Henceforth we refer to these as simply the close packed (cp) surface and the 100 surface, respectively. We work at the temperature k B T = 0.5 , where is the LJ potential well depth. As we are below the triple-point temperature k B T ≈ 0.65 [10], the liquid is less stable than the crystal, and so the crystal phase can nucleate on the surface. A crystal nucleus forming on the surface is shown in fig. 1.
The particles in our simulation can be divided into surface particles (S -coloured red in fig. 1), which are held in fixed positions throughout the simulation, and moving particles (M). Our results are all for 3520 moving particles. The surface consists of 3 layers of 20×22 particles, so that the system has 4840 particles in total. We define the lattice parameter of the surface a S as the nearest-neighbor distance between the surface particles (the red particles in fig. 1).
The interactions in our system are governed by a truncated and shifted LJ potential. We reduce the well depth of the cross-potential, MS / MM = 0.3, in order to inhibit rapid crystallisation of the moving particles. We choose σ MS such that a s = 2 1/6 σ MS , i.e. so that the minimum in the potential between moving and surface particles is equal to the lattice parameter of the surface. For simplicity, we define σ = σ MM and = MM , and work with σ and for the remainder of this letter.
In our simulations, the time evolution of the system is governed by NVT Monte Carlo (MC). To prepare the system, we initialise the fluid particles in random positions above the surface. We then "equilibrate" the system by evolving for 2×10 5 MC cycles, where each cycle consists of (on average) one attempted move per fluid particle. During this equilibration, the system is never observed to crystallise. Instead, the equilibration results in a liquid layer in contact with the surface, and in coexistence with the vapour. The liquid layer consists of 8-12 layers, which is significantly thicker than the crystal nuclei we study. We use the final particle configuration from the equilibration stage as the starting point for computing the nucleation rate.
In order to compute nucleation rates we use Forward Flux Sampling (FFS), a rare-event method developed by Allen et al. [15,16]. The order parameter in our FFS simulations is the size of the largest crystalline cluster, N cl , as identified by the local q 6 bond-order parameter introduced by ten Wolde and co-workers [12] and used in previous studies of LJ crystal nucleation [10,11,[17][18][19]. We refer to this cluster as the "nucleus" throughout this paper; thus the nucleus contains N cl particles.
Results and analysis. -When one crystal nucleates in contact with another, it has long been appreciated that the nucleation rate should be very sensitive to the match between the lattice of the surface and that of the nucleating crystal [2]. To study this we calculated nucleation rates on the cp and 100 surfaces as a function of the size of the mismatch between the lattice constant in the bulk nucleating phase, a B , and that of the surface, a S . This mismatch is quantified by [2] From MC simulations at constant pressure (NPT ) at k B T = 0.5 we computed a B ≈ 1.13σ, which is accurate to within 0.5%. We vary a S and hence δ, and plot the rate as a function of δ in fig. 2. We have also also performed simulations at k B T = 0.47 , for which the conclusions are the same. From fig. 2, we can see that the logarithm of the nucleation rate varies between ≈ −25 and −45 in our units. We can use this to estimate the nucleation barrier ΔG * , which can be written as ΔG * = −k B T log(k F F S /J 0 ), where J 0 is the kinetic prefactor that appears in the Classical Nucleation Theory (CNT) description [10]. If we assume that J 0 is approximately unity [1], the nucleation barrier is then ≈ 25-45k B T . Vonnegut and Turnbull [2] predicted that the nucleation rate is highest for δ = 0. This is not quite correct. As we can see in fig. 2, we find a maximum in the nucleation rate at a positive mismatch of ≈ 4% for the cp surface and ≈ 1% for the 100 surface.
To see if there is a perfect match between the surface and the nucleus at δ = 0, we calculated the average lattice parameter of the nucleus, x. This allows us to calculate the strain in the nucleus, i.e., the deviation from the bulk lattice constant: ξ = 100(x − a B )/a B . We obtained a value for x, and hence ξ, for each of the first five layers of the nuclei by averaging the nearest-neighbor separation for all pairs of neighbors (particles within 1.35σ) in that layer.
Our results for ξ are shown in fig. 3. We see that ξ is never zero, i.e., the spacing between the particles in the nucleus is never that of the bulk crystal. Our small nuclei have large surface-area-to-volume ratios and at surfaces (both crystal/liquid and crystal/crystal) the optimal value of the lattice constant will not in general be the bulk value. Isolated LJ clusters of comparable size have a lattice parameter that is strained ≈ 1% from the bulk lattice parameter [20].
Having determined that the interparticle spacing in the nucleus x is not the bulk lattice constant, we compare x with the surface lattice constant, a S . The two are equal when ξ = δ, so we plot the line ξ = δ in fig. 3. Where this line crosses the calculated strain in the first layer, the match between this layer and the surface is perfect.
For the cp surface, the ξ = δ line crosses the calculated strain at δ = 3, and this corresponds to the maximum in the nucleation rate. Therefore, for the cp surface the maximum in the rate does occur when the average lattice parameter of the first layer of the nucleus equals that of the surface. It is just that this does not occur when the lattice constant of the surface is equal to the bulk lattice constant, at δ = 0. For the 100 surface, the ξ = δ line again crosses the calculated strain at δ = 3. However, the maximum in the nucleation rate is at δ ≈ 1. Therefore, for the 100 surface the rate is not quite a maximum when the strain in the first layer is equal to the lattice constant of the surface. The complete surface has approximately double the linear dimension shown here (four times the area), and thus there is no interaction between particles in the nucleus across the periodic boundaries. For δ = −1 and 8, nucleation is found to be epitaxial for both surfaces. For the cp surface with δ = −13 (a), nucleation is not epitaxial. In this case the cp overlayer (and hence nucleus) is found to be oriented at a wide range of angles with respect to the surface -shown is an angle of ≈ 25 degrees. For the 100 substrate with δ = −13 (d), nucleation is epitaxial but the overlayer is cp, not 100.
The main point to make is that the finite nucleus does not have the lattice constant of the bulk crystal, therefore one would not expect the maximum in the nucleation rate to occur when the surface has the lattice constant of the bulk crystal. This fact will need to be taken into account when selecting or designing substrates for heterogeneous nucleation.
Looking at the first layers of nuclei on surfaces with a close match between nucleus and substrate lattices, fig. 4(b) and (e), we see that the two lattices are coherent. By coherent, we mean that the particles of the nucleus are all in equivalent positions (here the holes) in the surface lattice. As a S varies, in order to maintain coherence, the first layer of the nucleus has to deform, and indeed around δ = 0 in fig. 3, we see that there is definite slope to ξ as a function of δ, for δ near 0. However, the slope is less than one: on average the spacing between the particles in the nucleus does not perfectly track the surface lattice. We believe this is due to the finite size of the nucleus. Nuclei are always finite and this allows them to relax in from the sides to reduce the strain in the nucleus. We suggest that this feature of a slope less than one should be general, as all nuclei are finite, and so we propose to define coherent nucleation as being nucleation where ξ in the contact layer varies approximately linearly with δ. Note that when a nucleus is coherent, it will always be epitaxial, where by epitaxial we mean that the orientation of the substrate lattice fixes the orientation of the nucleus. Returning to fig. 2, we note that on both surfaces the rate roughly plateaus at negative values of δ. This is not what we would expect if nucleation is always coherent, as then the larger the value of |δ|, the larger the strain and the lower the rate. We can understand what is happening if we look at figs. 4 and 5.
In fig. 5 we have plotted histograms of the orientation angle the lattices of crystal nuclei make with the surface lattice. The orientation angle is defined as the angle between the unit cell of the surface and the first layer of the nucleus. We calculate this for a given nucleus by taking the dot product between the surface unit cell vectors and each vector r ij between two nearest-neighbour particles i and j. Figure 5(b) shows the distribution of nuclei orientations on the cp surface with δ = −1. It is very narrow: nucleation is epitaxial. Contrast this with the very broad distribution of angles for nuclei on surfaces with δ = −13 in fig. 5(a). For these large and negative mismatches the nucleation is no longer epitaxial. Here the nucleus is also not coherent, and so does not strain to fit the surface lattice, see also fig. 3.
The behavior on the 100 surface is different from that on the cp surface. For δ close to zero, the nucleus forms as expected with a 100 plane of the nucleus in contact with the 100 surface. However, as δ becomes more negative, the nucleus switches to forming with a cp layer in contact with the 100 surface. This is clear in fig. 4(d). Note that the nucleation is still epitaxial here. The close packed layer has a fixed orientation with a row of particles aligned with the surface rows. There is an intermediate "transition" region, −9 < δ < −4, where the first layer of the nuclei shows regions of sixfold symmetry and regions of fourfold symmetry. As shown in fig. 2, in this region the nucleation rate shows some evidence of a (shallow) minimum (at δ ≈ −6).
To get a better understanding of when nucleation is epitaxial, and which crystal plane of the nucleus forms in contact with the surface, we turn to the approach of Hillier and Ward [21] (henceforth HW). The approach of HW is based on the fact that the potential energy between a surface and an infinite crystalline overlayer is translationally invariant, and hence only depends on the orientation angle between the unit cell of the surface and that of the overlayer [22]. HW calculate a normalised interaction between the two crystal lattices, as a function of angle, and for different lattice planes, to determine both the favored lattice plane of the crystal to form in contact with the surface, and the favored angle. Here we go a little beyond this, and use the fact that in simulation we know the potential. So here we calculate the energy of interaction of an 8 × 8 particles overlayer (roughly the size of a nucleus produced from our FFS simulations) in contact with the surface. We fixed the lattice constant of the overlayer at the bulk value a B , although varying this does not affect our conclusions.
We find that this method correctly predicts the orientation of the nuclei that form for both surfaces and all mismatches δ we have studied. We illustrate this for a few particular cases in fig. 6. The top plot in fig. 6 shows the energy per particle of the overlayer, for a cp overlayer on a cp surface. Where we observe coherent nucleation (δ = 4), we see that there are deep energy minima at the orientations at which we see nuclei. However, for δ = −13 where nucleation is not epitaxial, we observe that the energy is a much flatter function of angle, with no deep minima. Thus, the HW approach appears to work here, in the sense that the presence or absence of deep minima is predictive of whether the nucleation is epitaxial or not 1 .
In the bottom panel of fig. 6, we have plotted the energy per particle for both cp (red curve) and 100 (green) layers in contact with a 100 surface with δ = −13. Again the HW result -predicting an epitaxial cp overlayer-agrees with our simulation findings, cf. fig. 4(d). In summary then, these simple potential-energy calculations successfully predict the loss of epitaxy on the cp surface, and the switch from nucleation with a 100 layer to a cp layer in contact with the 100 substrate. This success is encouraging. It suggests that the simple approach of HW could in some cases be used to predict epitaxy without resorting to numerically intensive computer simulations of nucleation.
Conclusion. -In conclusion, we used exact computer simulation to study the nucleation of one crystal on another. We found that, as expected from prior work, the nucleation rate is largest when the fit between the lattices of the nucleus and the surface is best. But we found that this was not when the lattice constant of the surface equalled that of the nucleating crystal in the bulk. The nucleus is a finite crystal with large surface-to-volume ratio, and so has a different lattice constant to the bulk crystal. These findings are presumably general, and will need to be borne in mind when selecting a substrate to induce nucleation of a crystal. A perfect match of the two bulk lattices will not necessarily maximise the nucleation efficiency. We also find that the simple method of Hillier and Ward [21] predicts well the orientation of the nucleus on the substrate. * * * We acknowledge EPSRC for funding (EP/J006106/1).