A molecular dynamics study of heterogeneous nucleation in generic liquid/substrate systems with positive lattice misfit

Nucleation plays a critical role in many natural and technological processes, and nucleation control requires detailed understanding of nucleation process at atomic level. In this study, we investigate the atomistic mechanism of heterogeneous nucleation in generic systems of liquid/substrate with positive lattice misfit (the solid has larger atomic spacing than the substrate) using molecular dynamics (MD) simulations. We found that heterogeneous nucleation process in such systems can be best described by a 3-layer nucleation mechanism: formation of the completely ordered first layer with an epitaxial relationship with the top surface of the substrate; formation of vacancies in the second layer to accommodate lattice misfit; and creation of a nearly perfect crystal plane of the solid in the third layer that demarcates the end of nucleation and the start of crystal growth. This 3-layer nucleation process creates a 2D nucleus (a plane of the solid phase), which contrasts with the hemisphere of the solid (a 3D nucleus) in the classical nucleation theory (CNT). It is expected that this 3-layer nucleation mechanism will provide new insight for nucleation control through effective manipulation of the liquid/substrate interface.


Introduction
Nucleation is the initial stage of any first order phase transformation and is a phenomenon widely spread in both natural and industrial processes [1,2]. Classical nucleation theory (CNT) [3][4][5] is commonly used to describe both homogeneous and heterogeneous nucleation processes despite a number of concerns about its validity [6][7][8]. While homogeneous nucleation is rarely observed in reality due to the presence of impurities in the liquid [1], so far we have little knowledge about the heterogeneous nucleation process at atomic level due to difficulties encountered by both computer simulations or direct experimental observations [9][10][11].
The recent epitaxial nucleation model [12] represents an alternative atomistic mechanism for the heterogeneous nucleation. It proposes that heterogeneous nucleation proceeds layer-by-layer through structural templating at the liquid/substrate interface. This has been supported by the latest results from both computer simulations [13][14][15][16][17][18] and experimental observations [19][20][21][22][23][24][25][26][27]. It has been reported that the liquid adjacent to the substrate can have significant atomic ordering even above the liquidus [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. The liquid may become layered within a few atomic layers at the interface, where the in-plane atomic ordering may exist [14,15]. With high resolution transmission electron microscopy (HRTEM), Oh et al found that the liquid Al exhibits the pronounced in-plane atomic ordering within 3 atomic layers adjacent to the α-Al 2 O 3 , where the atoms take hcp sequence of the substrate, not fcc sequence of the bulk solid Al [23,24]. The classical molecular dynamics (MD) simulation reveals that the 'prefreezing' crystalline Pb of approximately 2-3 atomic layers in thickness forms at the interface in the liquid Pb adjacent to (111) surface of solid Cu at 625 K [16]. In a recent MD study [18], we identified a 2-dimensional (2D) ordered structure at the generic liquid/substrate interface, using the local bondorder analysis [29], and we refer this phenomenon as prenucleation. Prenucleation is attributed to a structural templating mechanism, where the solid atoms in the underneath layer provide low energy positions for the solid Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. atoms in the new layer [18]. All these observations are consistent with the hypernucleation hypothesis [30] and the adsorption model [31].
The importance of structural templating can be further demonstrated by the facts that prenucleation facilitates heterogeneous nucleation while the 2D ordered structure is compatible with the new solid phase, or impedes heterogeneous nucleation while its structure is 'wrong'. For example, segregated Au atoms at the interface between AuSi eutectic liquid and Si (111) substrate exhibit a pentagonal atomic arrangement, which acts as the main barrier of the heterogeneous nucleation [27]. On the other hand, a (112) Al 3 Ti 2-dimensional compound (2DC) on the (0001) TiB 2 surface formed during grain refiner production makes TiB 2 particles extremely potent for the nucleation of α-Al, by decreasing the lattice misfit from f=−4.2% to 0.09% [32]. However, the addition of a small amount of Zr in Al melts results in the dissolution of the Al 3 Ti 2DC and the formation of Ti 2 Zr 2DC on the TiB 2 surface, which has a large misfit of −4.2% with α-Al, rendering TiB 2 particles impotent [33]. In addition, the experimentally measured undercooling for heterogeneous nucleation reveals that the nucleation potency of a substrate is highly dependent on the lattice misfit between the substrate and the solid [34,35].
Accommodating lattice misfit at the interface plays a critical role in the nucleation process and becomes an important step of heterogeneous nucleation [12]. The ab initio MD simulations reveals that the Al atoms exhibit a fcc-like ordering within 3 atomic layers at the interface between the liquid Al and Ti-terminated TiB 2 substrate at an undercooling of 2 K, and the growth of α-Al is frustrated due to the lattice misfit between solid Al and TiB 2 [36]. Recently, we found that in the generic system with negative lattice misfit (the solid has smaller atomic spacing than the substrate) heterogeneous nucleation completes deterministically within three atomic layers by structural templating to create a 2D template from which the new phase can grow [37]. In this 3-layer nucleation mechanism, a partial edge dislocation network forms in the 1st layer to accommodate largely the lattice misfit; the 2nd layer twists an angle via the formation of a partial screw dislocation network to reduce lattice distortion; and the creation of a crystal plane of the solid in the 3rd layer to template further growth. This 3-layer nucleation mechanism delivers a 2D nucleus, which contrasts with the cap formation process (3D nucleus) delivered through stochastic structural fluctuation in the liquid. However, a logical question would be whether this 3-layer nucleation mechanism is applicable to the generic systems with positive lattice misfit. The stress status in the 2D ordered structure at the interface is tensile and compressive, respectively, for the systems with negative and positive lattice misfit, and this difference could produce significant effects on the heterogeneous nucleation process [16,38].
The objective of this study is to study both prenucleation and nucleation processes in the generic systems with positive lattice misfit, using MD simulations.

Simulation approach
Similar to the case of the generic liquid/substrate systems with negative lattice misfit [37], here we have set up a generic liquid/substrate simulation system with positive misfit. The generic liquid consists of liquid Al, and generic substrates have a fcc lattice with a 〈111〉 surface orientation. We use Al atoms to construct the generic fcc substrates, which has an initial lattice parameter (a=4.126 Å) corresponding to the value obtained at the calculated melting point. The substrate atoms were fixed to a perfect fcc lattice configuration with a specified lattice parameter, a′, to achieve the desirable lattice misfit, f (0f8%). The definition of the lattice misfit can be found in Ref. [12]. This generic system has a number of advantages: (1) elimination of the effect of chemical interaction on heterogeneous nucleation; (2) systematic variation in lattice misfit; and (3) the substrates can mimic dissimilar materials with much higher melting temperatures (T l ) (e.g., TiB 2 with T l =3498 K [39]). More details of such generic simulation system are given in Ref. [37].
The dimension of the substrate is 96[¯] 112 ×66[¯] 110 ×5[ ] 111 and the initial dimensions of the liquid is optimized to minimize the strain between the substrate and the liquid, with total numbers of atoms between 23550-26112. A vacuum region is inserted with periodic boundary conditions in the z ([111])direction, and the extent of the vacuum region is 60 Å. Periodic boundary conditions are also imposed in the x ([¯] 112 ) and y ([¯] 110 ) directions. It has been confirmed that simulations will produce only minor differences in the degree of atomic ordering at the interface for the systems with relaxed and fixed substrates [40].
The EAM potential for Al, developed by Zope and Mishin [41], was used to model the interatomic interactions. The predicted melting temperature is 870±4 K for pure Al with the EAM potential [41]. MD simulations are performed using the DL_POLY_4.08 MD package [42]. The equations of motion are integrated by means of the Verlet algorithm with a time step of 0.001 ps, and the Berendsen NVT ensemble is used for the temperature control. The liquid sample is prepared by heating the system to a temperature of 1400 K with steps of 50 K, each lasting 100,000 MD steps.
The equilibrated configuration of the liquid was cooled from 1400 K to the desired temperature with steps of 50 K, 5 K and then 1 K, and at each temperature step the system was allowed to run for 1,000,000 MD steps to equilibrate. The nucleation temperature, T n , was determined by monitoring variation in total energy and trajectory of the system during the equilibration. Firstly, we determine the temperature, T 1 , at which nucleation occurred, with a temperature step of 50 K. The configuration of the system equilibrated at T 1 +50 K is used as initial configuration for the next simulation with a temperature step of 5 K to determine the next temperature, T 2 , where the nucleation occurred. Similarly, the configuration of the system equilibrated at T 2 +5 K with a temperature step of 5 K is used as initial configuration for the next simulation with a temperature step of 1 K, and this approach allows us to pinpoint T n within ±1 K.
Atomic layering at the interface is quantified by the atomic density profile, ρ(z), where z is the distance along the z ([111]) direction of the simulation system [14]. The atomic ordering in each layer of the liquid at the interface is characterized by the z-dependent in-plane order parameter, S(z) [43]. The details for the calculations of ρ(z) and S(z) can be found in Ref. [18].

Prenucleation
Following our previous work on prenucleation in the generic liquid/substrate system with negative misfit [18,37], we have investigated the prenucleation phenomenon in the liquid/substrate system with positive misfit. Here we use the system with 4% misfit as an example to demonstrate the prenucleation phenomenon in detail. Figure 1(a) shows a snapshot of the simulation system with f=4% equilibrated at T=1000 K, which is 130 K above the liquidus determined by MD calculation. The liquid adjacent to the interface becomes layered within about 6 atomic layers, as shown by the quantified atomic density profile in figure 1 shown the time-averaged atomic positions of the first layer (L1) on the top of those of the surface layer of the substrate (L0), the second layer (L2) and the third layer (L3) at the interface, respectively. Figure 1(c) suggests that L1 has a completely ordered atomic arrangement and is in an epitaxial relationship with the substrate surface (L0). L2 has a mixed structure with ordered and disordered regions, and the L3 is almost completely disordered. This suggests that there exists a 2D ordered structure in the liquid adjacent to the interface, which is in qualitatively agreement with that in the system with negative misfit [18,37].
We also have examined the effect of temperature and lattice misfit on the degree of prenucleation. Figures 2(a) and (b) show the quantified peak density, ρ p (z), and in-plane order parameter, S(z), in the system with f=8% as a function of temperature, respectively. With a nearly 500 K decrease in temperature, both ρ p (z) and S(z) increases only slightly, suggesting that temperature only have a moderate effect on atomic ordering in the liquid adjacent to the interface. Figures 2(c) and (d) represent the quantified ρ p (z) and S(z) in the liquid equilibrated at 1000 K as a function of lattice misfit, respectively. These results suggest that lattice misfit has no effect on atomic layering, but has a strong influence on in-plane atomic ordering. In-plane atomic ordering increases sharply with decreasing lattice misfit (figure 2(d)), particularly in the first 2 atomic layers. This is consistent with the results for the simulation system with negative misfit [18,37]. We can conclude that while temperature only have a moderate effect on prenucleation, lattice misfit has a strong influence on atomic ordering in the liquid adjacent to the interface regardless of the nature of the lattice misfit, positive or negative.

Nucleation
Now, we take the system with f=8% as an example to investigate the atomistic mechanisms of heterogeneous nucleation in the generic systems with positive lattice misfit. Heterogeneous nucleation was observed by monitoring the total energy, E t , of the simulation system as a function of the lapsed time, t, with increasing the undercooling from T l . Figure 3 displays E t as a function of rescaled time, t (taking the start of nucleation as t=0 ps) for the system with f=8% during the simulation at T=739 K. E t remains unchanged before t=0 ps, and drops dramatically from 0 ps to 180 ps, suggesting that the nucleation occurs at T n =739 K for this simulation system. a and b in figure 3 marks the onset and the end of the solidification process, respectively. Figure 3 suggests that at the nucleation temperature heterogeneous nucleation is a spontaneous down-hill process without the need to overcome any energy barrier.
To understand the atomistic mechanism of heterogeneous nucleation, we have examined the evolution of atomic ordering in the liquid adjacent to the interface at the nucleation temperature (T n ). Figure 4 shows the time-averaged atomic positions of the first 3 atomic layers (L1, L2 and L3) of the simulation system with f=8% from t=−100 ps to −10 ps during the simulation at T n =739 K. During this period, L1 remains completely ordered, and the L2 consists of the ordered and disordered regions, and the L3 is largely disordered, with only a few ordered regions being visible. It is noted that the size and position of the ordered regions in L3 are dynamic and change from time to time, suggesting that the simulation system is still in the stage of the prenucleation before t=0 ps. Figure 5 shows the evolution of the 2D ordered structure in the liquid adjacent to the interface during heterogeneous nucleation for the simulation system with f=8% at T n =739 K. For each time, we show the front view (a snapshot) and top views (time averaged atomic positions) of L2 and L3. The onset of nucleation (t=0 ps) is symbolized by the stabilization and subsequent growth of a particular ordered region in L2 and L3, as marked by the circle in figure 5(a3). Figure 5(a) (t=0 ps) represent the maximum degree of atomic ordering at the interface that can be achieved by prenucleation, and therefore is treated as a precursor for the heterogeneous nucleation process. Building on the precursor, this ordered 2D structure expands (in 2D) continuously with time (marked by the circles) and reaches 4.3 nm at t=40 ps ( figure 5(c3)), which matches the critical size of nuclei at T n =739 K determined by the classical nucleation theory (CNT) and marks the end of the nucleation. It should be pointed out that for t<40 ps, the growth is in 2D, i.e., lateral growth of the precursor on the planes parallel to the interface. After nucleation, further growth of the nuclei will be 3D and leads to the formation of a hemisphere at t=70 ps, as highlighted by the half circle in figure 5(d1). Beyond this point, growth enters spherical growth stage. In this paper we will confine our discussions to prenucleation and heterogeneous nucleation, and leave cap formation and spherical growth to other publications.
The atomic spacing, d a , of the ordered regions within an individual layer at the interface is calculated by averaging the distances between all the nearest neighboring atoms. Figure 6(a) shows the variation of d a of the individual layer at the interface as a function of time for the system with f=8% during the simulation at T n =739 K. d a of L1 is about 2.7 Å, which is very close to 2.68 Å of the substrate, and it is constant with the observation that the L1 is epitaxial to the surface layer of the substrate ( figure 4). d a is then increased to 2.83 Å in L2 and 2.92 Å in L3. Beyond L3 d a becomes constant, being the atomic spacing on (111) of aluminum. The   figure 6(a) suggest that L2 is a transitional layer between the highly compressed L1 and the stress-free L3 and beyond.
In order to understand the mechanism for accommodation of lattice misfit at the interface we have calculated f a , which is the lattice misfit accommodated by a given layer relative to the previous layer. f a as a function of the atomic layers away from the interface for the systems with f=8% at 739 K and t=1 ns is plotted in figure 6(b). f a is negligible between L1 and the surface layer of the substrate (L0), and about 8.5% between the L2 and L1, and approaches 0% beyond L2. Figure 6(b) indicates that all the apparent lattice misfit of the simulation system is accommodated by L2. Further analysis suggests that the accommodation of lattice misfit in L2 is attributed to the formation of vacancies in L2 during the nucleation process. Figure 7(a) shows the snapshots of L1, L2 and L3 for the systems with f=8% at 739 K at t=1 ns. Examination of atomic arrangements in figure 7(a) revealed a three-layer process for heterogenous nucleation: (1) L1 is a nearly perfect crystalline plane that is in epitaxial relationship with L0; (2) L2 contains 17.3% vacancies; and (3) L3 is a nearly perfect crystalline plane of the solid phase with some equilibrium vacancies.
We have investigated the applicability of this three-layer process to systems with different lattice misfit. Figure 7 shows the snapshots of L1, L2 and L3 for the systems with f=2%-8% at t=1 ns and at their corresponding nucleation temperatures. It is interesting to note that independent of lattice misfit of the system L1 is always in epitaxy with L0 and L3 is always a nearly perfect (111) plane of the solid containing a small amount of vacancies (assuming equilibrium vacancies). The only observable change with lattice misfit is the level of vacancies in L2, which increases with increasing apparent lattice misfit ( figure 7). This is confirmed by the quantified levels of vacancies in L2 as a function of the apparent lattice misfit ( figure 8). Nucleation undercooling, ΔT n , is calculated as ΔT n =T l − T n , and ΔT n as a function of f for all the systems under the same simulation conditions is plotted in figure 9. The data for the systems with negative misfit [39] are also included for comparison. ΔT n increases significantly with increasing apparent lattice misfit in the studied range of lattice misfit (0f8%), regardless of the sign of the misfit being negative or positive. If ΔT n is used as a quantitative measure of nucleation potency of a substrate (the smaller the ΔT n , the higher the nucleation potency), the nucleation potency decreases significantly with the increase of lattice misfit. This suggests that the lattice misfit is an important factor in determining the potency of the nucleating substrates. (b) during nucleation (t=20 ps), (c) completion of nucleation (t=40 ps); (d) during the cap formation (t=70 ps). The ordered 2D structure highlighted by a circle in L3 expands in 2D with time, reaches a disk diameter of 4.3 nm at t=40 ps (c3), which corresponds to the critical cluster diameter defined by CNT at ΔT n =131 K. This is followed by isothermal cap formation ((d1)-(d3)).

Discussion
In a recent study [37], we defined heterogeneous nucleation from an atomistic point of view. Heterogeneous nucleation process starts with a precursor (a 2D ordered structure) at the liquid/substrate interface created by prenucleation [18], proceeds layer-by-layer through structural templating and finishes by creating a crystal plane of the solid, from which further growth can be templated. In the same study [37], we demonstrated a three-layer mechanism for heterogeneous nucleation in a generic system with negative lattice misfit (the solid has smaller atomic spacing than the substrate): L1 reduces lattice misfit by forming a set of edge dislocations; L2 twists an angle to reduce the lattice distortion by forming a set of screw dislocations; and L3 is a perfect plane of the solid, which will template further growth of the solid. This three-layer nucleation process can be also described as a transformation from a liquid/substrate interface to two new interfaces: a substrate/solid interface and a liquid/ solid interface.
In this study we revealed an atomistic nucleation mechanism in the generic systems with positive lattice misfit (the solid has larger atomic spacing than the substrate), which is schematically illustrated in figure 10. At the prenucleation stage (T>T n ), the lattice of the substrate induces the prenucleation in the liquid at the interface by structural templating, resulting in a 2D ordered structure (figure 5(a)) in which the first layer (L1) is completely ordered with an epitaxial relationship with the surface plane of the substrate (figure 4). During heterogeneous nucleation at nucleation temperature, L2 is templated by L1 and accommodates all the lattice misfit by formation of vacancies randomly populated in L2 (figure 7), and this is followed by the formation of L3 that is a crystal plane of the solid. L3 templates the growth of further layers by advancing the liquid/solid interface into the liquid, and solidification enters growth stage. Compared with the nucleation mechanism in the systems with negative lattice misfit [37], the nucleation mechanism in the systems with positive misfit can be also described as a three-layer nucleation mechanism.
The major difference between the two systems is how the lattice misfit is accommodated: negative lattice misfit is accommodated by the formation of edge dislocations in L1 followed by screw dislocations in L2, while positive lattice misfit by epitaxy in L1 followed by the formation of vacancies in L2. Although further investigation is required to understand the exact cause of this difference between systems with positive and negative misfit, it is possible that the nature of the local stress distribution across the interface may play an important role in determining the mechanisms for misfit accommodation. In the system with a negative lattice misfit, the ordered atoms induced by the substrate is under tensile stress with a possible maximum stress being vertical to the 〈011〉directions. This suggests that it is easier to introduce extra rows of atoms along the 〈011〉directions, leading to the formation of 〈011〉edge dislocations. However, in the case of systems with positive lattice misfit, the ordered atom induced by the substrate are under compression stress along all the directions due to their epitaxy relationship with the substrate lattice. In this case, it would be easier to take randomly individual atoms out (forming vacancies) to release the tensile stresses.
This study suggests that the fraction of vacancies in the second layer of the solid increases with increasing apparent lattice misfit in the system (figure 8). Here we derive the relationship between the apparent lattice misfit ( f, or accommodated lattice misfit, f a ) and vacancy fraction ( f v ) in L2.  Insert the definition of vacancy fraction ( f v =(n sub -n m )/n sub ) into equation (2) we obtain the following relationship between f a and f v : Predictions by equation (3) agrees well with the results from MD simulations, as show in figure 8.
In addition, formation of vacancies as an mechanism to accommodate positive lattice misfit was reported in our previous ab initio study of the liquid Mg/MgO system with {111} MgO surface termination, where the ordered Mg atoms in the liquid induced by the {111} MgO substrate containing significant amount of vacancies [44].
Finally, we should point out that in this study nucleation is defined as a process to form a crystal plane of the solid (a 2D nucleus) which can template further crystal growth. This is obviously different from the formation of a spherical cap (a 3D nucleus) as defined by CNT. In CNT, the nucleation potency (ΔT n ) of a substrate is described as a function of contact angle (θ). θ, as a function of the relevant interfacial energies, is the outcome of a continuum analysis of heterogeneous nucleation process. Due to the lack of reliable data for interfacial energies, θ can not provide further insight for nucleation control. In contrast to CNT, the new definition of nucleation (2D nucleus) suggest that nucleation potency of a given substrate is closely related to lattice misfit [37], chemical  interaction between the liquid and the substrate [40,44] and the atomic level surface roughness [33,45]. This knowledge provides further insight for nucleation control, in particularly through manipulation of the liquid/ substrate interface via chemical segregation. For instance, TiB 2 has a hexagonal crystal structure and a hexagonal platelet morphology with the (0001) planes as its major surface termination. The misfit between (0001) TiB 2 surface and {111} Al is -4.2% at 660°C [32], suggesting that TiB 2 is not potent for the nucleation of Al. Segregation of Ti at the TiB 2 /liquid Al interface leads to the formation of a mono atomic layer of (112) Al 3 Ti 2DC, which reduces the |f| from 4.2% to 0.09% and increases significantly the nucleation potency of TiB 2 [32]. However, addition of 500ppm Zr into Al melt results in the formation of (0001) Ti 2 Zr 2DC on (0001) TiB 2 surface, which replaces the original (112) Al 3 Ti 2DC [33]. The (0001) Ti 2 Zr 2DC has not only a large misfit with Al but also an atomically rough surface, making such TiB 2 particles impotent for the nucleation of solid Al.

Summary
In this work, we used MD simulation approach to investigate the effect of positive lattice misfit on both prenucleation and heterogeneous nucleation processes. We found that similar to the systems with negative lattice misfit, prenucleation in the systems with positive lattice misfit is promoted by reducing both temperature and lattice misfit. More importantly, our study also revealed that heterogenous nucleation process is complete with 3 atomic layers (the 3-layer mechanism): complete ordering of the first layer in epitaxy with the substrate; formation of vacancies in the second layer to accommodate lattice misfit; and formation of a nearly perfect crystal plane of the solid (the 2D nucleus) to complete nucleation. This 3-layer nucleation mechanism for the systems with positive misfit is qualitatively similar to that for the systems with negative misfit, but differ in the detailed mechanisms for accommodating lattice misfit; the former accommodates misfit through formation of Figure 10. Schematic illustration of the 3-layer atomistic mechanism for heterogeneous nucleation in the systems with positive lattice misfit. The 1st layer (L1) formed during prenucleation is completely ordered and is epitaxial to the surface layer of the substrate; at the nucleation temperature, the epitaxial L1 templates the 2nd ordered layer (L2) which contain certain amount of vacancies to accommodate lattice misfit; and L2 templates a nearly perfect crystalline plane of the solid in L3, the 2D nucleus, to template further growth.
vacancies while the latter through formation of dislocations. Such detailed atomistic mechanisms may provide insight for nucleation control through manipulation of the liquid/substrate interface.