C

rystallization from the melt (Fig. 1) is a pervasive process in industry, from metal casting for structural applications to the Czochralski process for semiconductor wafer growth for electronics. It is important to control and understand the crystal growth process because it is at this stage that the material's microstructure morphology is created, which in turn defines the material's properties. Consequently, a great deal of effort has been put into understanding the complex interplay between structure, thermodynamics, and kinetics that governs the process of crystal growth [1][2][3] . This has led to a mechanism-based understanding of crystallization [4][5][6] in terms of the microstructural elements of the crystallite being formed. For example, the character of the solid surface in contact with the liquid is known to affect the growth rate, with atomically rough surfaces leading to faster growth rates than flat low-index surfaces and their vicinals. Considerably less attention has been put in understanding the effects that the liquid adjacent to the solid-liquid interface has on the process of crystal growth.

Atomic events leading to crystal growth are thermally activated processes taking place in the free-energy landscape illustrated in Fig. 2a. The rate of crystallization is proportional to expðÀβΔE a Þ, while the melting rate is proportional to exp½ÀβðΔE a þ ΔμÞ, where ΔE a is the activation energy for solidification, Δμ is the difference in chemical potential between the liquid and solid phases, β -1 ≡ k B T, and k B is the Boltzmann constant. The balance of these two rates results in the following equation for the overall growth rate:

where kðTÞ k 0 expðÀβΔE a Þ is known as the kinetic factor. In this model, known as the Wilson-Frenkel 7,8 (WF) model, the activation energy for solidification is taken as the energy barrier for diffusion in the liquid, ΔE a = ΔE d , because crystallizing atoms must undergo the same self-diffusion process that occurs in the associated liquid phase. It is often found that the WF method cannot quantitatively predict results from simulations or experiments 9,10 . This notorious discrepancy, while largely unsolved, has been attributed to changes in mobility of the supercooled liquid in the vicinity of the crystal interface that would cause ΔE a > ΔE d , but no physical mechanism has been demonstrated to explain the origin of this effect.

Here we employ atomistic simulations and machine learning (ML) together to show that the solid-liquid interface induces partial ordering of the nearby liquid during crystal growth. Our approach is successfully applied to two different families of materials: semiconductors and metals. We find that the interfaceinduced ordering (IIO) affects the mobility of liquid atoms and thus slows down the crystal growth kinetics. The physical mechanism behind the IIO is explained and we demonstrate that by accounting for this effect it is possible to derive predictive models for crystal growth.

Results

Crystal growth simulations. We performed molecular dynamics (MD) simulations of crystalline silicon growth from its melt employing a simulation geometry akin to laboratory experiments of crystal growth: a crystalline seed is introduced in the liquid and its growth is monitored over the course of the simulation (see Fig. 1 and Supplementary Video 1). This setup allows the different microstructural elements of the growing crystallite to interact naturally (see Supplementary Videos 2 and 3), as they would in a crystal growth experiment. For this geometry Δμ = ΔG-κγ/ρ s , where ΔG is the difference in free energy between the liquid and the crystal, and the second term is due to the Gibbs-Thomson effect, with ρ s being the density of the solid, γ the interfacial free energy, and κ ¼ 2=R eff is a geometrical factor where R eff is the effective crystallite radius. All the above parameters of the WF model were computed (Fig. 2b and c) in order to compare the model predictions against simulation results (for calculation details see "Methods" section and Supplementary Notes 5 and 6). The comparison between model and simulations is shown in Fig. 3a, where it is evident that the WF model does not predict the growth rate for temperatures it was not fitted to.

Machine learning encoding of crystallization events. Historically, simpler simulation geometries have been favored as a way to isolate certain microstructural elements, which are then probed separately [10][11][12] . Our use of the geometry shown in Fig. 1 makes the simulation more physically relevant at the expense of greatly diminishing the amount of information that can be inferred due to the lack of a crystal growth model that accounts for all microstructural elements present and their respective interactions. Moreover, it also becomes challenging to decipher the atomic events at play due to the sheer complexity of the environment that atoms are embedded in. Here these obstacles are overcome by employing ML algorithms to systematically encode and classify the structure surrounding liquid atoms during crystallization events. Our approach builds on recently proposed ML strategies for the construction of a structural quantity (namely softness S) that captures the propensity for atomic rearrangements to occur in disordered atomic environments, such as in glasses 13,14 and inside grain boundaries 15 .

The structural characterization of local atomic environments is realized by assigning to each atom i a local-structure fingerprint x i constructed from a set of 21 radial structure functions 13,16 GðrÞ, as illustrated in Fig. 4a. Furthermore, atoms are labeled into three possible categories according to their first-neighbor's arrangement: liquid and crystal atoms have arrangement patterns statistically identical to the bulk liquid and bulk crystal, respectively. Meanwhile, crystallizing atoms have arrangement patterns intermediary between the other two labels (see "Methods" section and Supplementary Note 2 for more details). t = 0 ps t = 981 ps t = 1659 ps t = 2090 ps Fig. 1 Crystal growth simulations. Snapshots of a crystal growth simulation of silicon using molecular dynamics. The system initially contains a small crystalline seed (yellow atoms) surrounded by liquid (transparent blue atoms). Shown in the last frame is the dislocation network 46,62,63 formed during the growth process (edge dislocations are colored blue while screw dislocations are shown in red). See Supplementary Video 1 for the complete video and Supplementary Video 3 for dislocation network evolution.

It is possible to observe how these three groups of atoms are spread in the R 21 -space of local-structure fingerprints x i with the help of an algorithm known as principal component analysis (PCA). With this method, a dimensionality reduction transformation is performed to create a R 2 representation of the R 21 data, as shown in Fig. 4b. Superimposed in this figure is also the trajectory of an atom that undergoes crystallization over the course of the simulation. Atoms assume varied local-structure fingerprints x i depending on both the surrounding liquid structure and the nearby interface morphology. In order to quantify these variations in microstructure we proceed as follows. First, an ML algorithm known as support vector machine [17][18][19] is employed to find the hyperplane that optimally separates the crystallizing atoms from the liquid atoms in the R 21 -space of x i . Then, the distance of each atom i from the hyperplane (S i , known as softness [13][14][15] ) is measured: atoms with S i >0 lie on the crystallizing side of the hyperplane, while S i <0 atoms lie on the liquid side. This approach is found to correctly classify liquid and crystallizing atoms with an accuracy of 96%. It is important to realize that S is not an order parameter because it was not designed to track the change from the liquid to the solid phase. Instead, S measures the propensity of an atom in the liquid phase to undergo the process of crystallization.

Shown in Fig. 4c is a simulation snapshot with atoms colored according to their softness value (see also Supplementary Videos 4 and 5). In this figure S is seen to capture the structural signs of dynamical heterogeneity in the supercooled liquid far from the crystal, with clear indications of strong spatial correlations. These fluctuating heterogeneities have recently been shown to be preferential sites for crystal nucleation 20,21 . Thus, Fig. 4c establishes that S is indeed capable of capturing subtle signs of structural ordering in liquids.

Local-structure dependent (LSD) crystal growth model. It is possible now to decompose the crystal growth rate of Fig. 3a as a function of the local structure using S. For example, in Fig. 3b it a b 0.5 0.4 0.3 0.2 0.1 0.0 1100 1200 1300 Molecular dynamics simulations Data used for model parametrization Wilson-Frenkel model Local-structure dependent model Pulsed-laser melting experiments (Ref. 64) 1400 Interface temperature [K] Growth rate [atoms per nm 2 ps] Growth rate [atoms per nm 2 ps ] 1500 1600 -1.5 10 -1 10 -2 10 -3 10 -4 -1.0 -0.5 0.0 0.5 1.0 1.5

Fig. 3 Local-structure dependent crystal growth model. a Crystal growth rate of silicon versus interface temperature. Both models (Wilson-Frenkel and local-structure dependent model) were parametrized using only the data from simulations at T ≥ 1388K. Notice that the growth rate is given per unit of effective area of the total crystallite (see Supplementary Note 4) and the error bar represents the 95% confidence interval around the mean (see "Methods" section). Data for laser-pulsed melting experiments was extracted from Galvin et al. 64 . b Growth rate at T ≈ 1233K decomposed as a function of the local structure, as encoded by S, surrounding crystallizing atoms. Notice that the growth rate varies almost four orders of magnitude as S changes. is shown that the total growth rate at T ≈ 1233 K varies by almost four orders of magnitude as the local structure changes. For this reason, we propose to address the limitations of the WF model by taking into account the local structure surrounding the crystallizing atoms through an explicit dependence on S: rðT; SÞ ¼ kðT; SÞ 1 À exp ÀβΔμðT; SÞ ½ f g ;

with ΔμðT; SÞ ¼ ΔGðTÞ À κγðSÞ=ρ s . Indeed, accounting for the information about the local structure contained in plots such as Fig. 3b results in a crystal growth model with predictive capabilities, as shown in Fig. 3a (see Supplementary Notes 3 and 4 for details on the model calculation). Notice how the LSD model is capable of predicting the growth rate for a wide range of temperatures (i.e. T < 1388 K) not included in the model parametrization. In particular, the experimentally measured growth rate and its slope show much better agreement with our LSD model, Eq. ( 2), than with the WF model, Eq. (1). In the Supplementary Note 3 we show that the variables introduced by the dependence on S are not independent parameters. Thus, the improved reproduction and prediction of simulation results cannot be attributed simply to Eq. ( 2) exhibiting higher capacity or flexibility in modeling complex relationships when compared to Eq. ( 1). We now turn to investigate the ramifications of the LSD model rðT; SÞ and uncover the source of its predictive capabilities. The kinetic factor kðT; SÞ, shown in Fig. 5a, is observed to be a strong function of the local structure, varying by as much as three orders of magnitude with S. For each value of S the kinetic factor shows an Arrhenius-like temperature dependence kðT; SÞ ¼ k 0 ðSÞ exp½ÀβΔE a ðSÞ:

This striking outcome suggests that each value of S corresponds to a thermally activated and independent crystallization channel with well-characterized energy scale. Such a picture is reminiscent of our traditional understanding of crystallization in terms of the solid-liquid interface morphology, with different values of S encoding the influence of different microstructural elements. But here S encodes more than just the crystal local microstructure: it also encodes the variation in the structure of the liquid. The variation of liquid properties with local structure is reflected in the dependence of the activation energy barrier of these crystallization channels with S shown in Fig. 5b, where it can be seen that ΔE a ðSÞ varies over 1 eV with S. Additionally, the activation energy decreases monotonically with S and seems to approach the energy barrier for diffusion ΔE d (Fig. 2b) asymptotically. Hence, the mobility of liquid atoms close to the solid-liquid interface seems to vary greatly, from a negligible change (ΔE a ≈ ΔE d ) compared to bulk liquid to a dramatic reduction in mobility due to the increase in ΔE a . This change in the liquid structure due to the presence of the solid-liquid interface is known as IIO, Supplementary Fig. 5b shows that the structural change does indeed lead to local ordering.

Varying S also has a pronounced effect on the Arrhenius prefactor k 0 ðSÞ, as indicated in Fig. 5c, which decreases by three orders of magnitude with S. Because ln ½k 0 ðSÞ can be interpreted as the product of an entropic contribution 15 to the free energy barrier and a term involving the population of crystallizing atoms with softness S, the observed decrease in the prefactor indicates that there are less rearrangement pathways leading liquid atoms to the activated state (i.e. to crystallization) as S increases. Hence, Fig. 5b and c together indicate that from all observed localstructure arrangements surrounding crystallizing atoms, only very few lead to low-energy barriers. Additionally, Fig. 3b indicates that these few channels with low-energy barriers are the ones contributing the most to the overall growth rate.

Next, we examine how the free energy of the solid-liquid interface to which atoms attach varies with S, which should give us a glimpse of the microstructure at the crystallite surface. Figure 6a shows that γðSÞ decreases monotonically with S, starting at large values-corresponding to high-index interfacesand reaching interfacial free energy values characteristic of lowindex interfaces in silicon. This finding implies that the decrease in Arrhenius prefactor k 0 ðSÞ with softness (Fig. 5c) leads to fewer rearrangement pathways because crystallization events with large positive values of S happen at low-index surfaces and their vicinals, which naturally offer less crystallization sites than highindex interfaces. Despite the scarcity of crystallization sites offered by low-index interfaces and their vicinals, Fig. 3b shows

a b c 4 2 0 -2 -4 -6 -8 -5 0 Second principal component First principal component Liquid Crystal Crystallizing Atom trajectory 5 -1.5 1.5 r 3 r 4 r 2 r 1 r cut that they contribute the most to the overall growth rate, with 70% of all atoms attaching to interfaces with S ≥ 0:75. This observation is confirmed by direct measurement of the distribution of crystal surfaces to which crystallizing atoms attach: Fig. 6b reveals strong preferential attachment to a wide variety of (111) vicinals. The high-intensity spot around (435) corresponds to step-step separation distances from 15 to 24 Å (Fig. 6c), indicating that the majority of the crystallization events take place on vicinal surfaces with well-separated steps, which is exactly what is expected for silicon 22 . Figure 6b also shows a smaller amount of events occurring at high-index interfaces, further validating the above observations. Notice in Supplementary Video 1 that the crystallite also exhibits signs of rough interfaces, thus it is possible that a small fraction of the identified high-index are actually rough.

Applicability to a different family of materials. In order to verify that our approach in creating LSD predictive models of crystal growth is not particular to silicon (or semiconductors) we apply it in the development of a crystal growth model for an elemental metal, namely copper (see "Methods" section and Supplementary Note 7 for simulation details). The resulting model is shown in Fig. 7a, where it can be seen that the LSD model of copper also correctly predicts the growth rate at temperatures at which it was not parametrized on (i.e. it is a predictive model), while the WF model is not capable of reproducing simulation results at temperatures to which it was not fitted, similarly to what was observed for silicon in Fig. 3a. Moreover, analysis of the parameters of the LSD model of copper (Fig. 7b) shows that all parameters present the same trend with S as observed for silicon, including the Arrhenius behavior for the kinetic factor (shown in Supplementary Fig. 19c). The major difference observed between the LSD models of silicon and copper is that k 0 ðSÞ and ΔE a ðSÞ assume much larger values for copper. We attribute this to the predominance of rough interfaces in metallic systems. In contrast to semiconductors, interfaces in metallic systems typically do not advance by the lateral motion of steps. Instead, metal interfaces often advance by atomic attachment directly on top of them, leading to growth normal to the interface itself. This growth mechanism is reflected in Fig. 7a (inset), where it is seen that atomic attachments occur directly on (001) and (111) interfaces-leading to normal growth -instead of vicinals of (111) as observed for silicon in Fig. 6b (compare also Fig. 1 to Supplementary Fig. 14 and Supplementary Video 1 to Supplementary Video 7). Normal growth leads to the formation of atomically rough interfaces that offer a much larger amount of atomic disorder than well-structured high-index interfaces. Hence, rough interfaces present a larger availability of sites for liquid atoms to attach, leading to much higher values for k 0 ðSÞ due to the numerous atomic pathways leading to crystallization. The predominance of rough interfaces also explains the larger values of ΔE a ðSÞ observed for copper, but this explanation will be postponed until the "Discussion" section, where the connection will be discussed in the light of the effects of IIO on crystal growth.

Discussion

Solid-liquid interfaces in equilibrium are known to affect the structure of the nearby liquid by imparting some amount of order on it [23][24][25][26][27][28][29][30][31] . Here, we have established that the IIO of the liquid also occurs during the process of crystal growth-a dynamic situation in which the solid-liquid interface is not in equilibrium. The observed IIO seems to decrease the mobility of liquid atoms through changes in the activation barrier for crystallization ΔE a (Fig. 5b), effectively slowing down the crystallization kinetics.

Comparison of Figs. 5b and 6a reveals that the IIO of the liquid is anisotropic, i.e. it depends on the interface orientation and microscopic levels of roughness. The trend (illustrated in Fig. 8) is such that low-index surfaces and their vicinals (corresponding to large positive S in Fig. 6a) cause weak ordering, resulting in smaller activation energies (i.e. ΔE a ðSÞ close to the energy barrier for atomic diffusion in the liquid bulk ΔE d ). Meanwhile, highindex interfaces (negative S in Fig. 6a) cause strong ordering of the liquid, which becomes rigid and results in activation energies much larger than ΔE d . However, even in the case of strong IIO the activation energy (1.75 eV) is still much smaller than the ≈4.6 eV 32 barrier for vacancy-mediated self-diffusion in crystalline silicon. This indicates that the structural order of the liquid affected by IIO is nowhere as substantial as crystalline order.

The physical cause of the IIO anisotropy is that the interaction between the crystal surface and the liquid is mediated by the amount of dangling bonds on the crystal surface. Thus, strong liquid ordering (and slower mobility) is observed at high-index interfaces because these interfaces interact more strongly with the liquid, since they present more dangling bonds when compared to low-index interfaces and its vicinals. This mechanism is illustrated in Fig. 8a and b, while its effect on the free-energy landscape of the system is illustrated schematically in Fig. 8c. Notice that this mechanism also explains why copper has much larger values of ΔE a ðSÞ (Fig. 7b) while having similar energy barrier for diffusion in the liquid ΔE d : rough interfaces are predominant in copper and these interfaces have much stronger interactions with the liquid when compared to low-index and vicinal surfaces, which are predominant in silicon.

Dynamical heterogeneities present in the liquid (Fig. 4c) also affect the coordination of atoms 20,33,34 . For this reason, it is reasonable to expect that they contribute to the ≈1 eV dispersion in ΔE a ðSÞ observed in Fig. 5b. Nonetheless, there is no evident reason to believe that their effect is anisotropic since dynamical heterogeneities have origin in random thermal fluctuations.

In conclusion, we have discovered that the IIO of liquids strongly affects the process of crystal growth in metals and semiconductors. It is found that the modified structure of the liquid nearby solid-liquid interfaces reduces the mobility of liquid atoms, an effect shown to be essential in order to build a predictive model of the growth rate temperature dependence. Indeed, the construction of such predictive model was only possible by identifying and incorporating in the model the family of all thermally activated events-each with its own energy scaleleading liquid atoms to the crystal phase. Our work elevates the

Low-index or vicinal surface High-index or rough surface 1.2 a b c ( ) (100) (101) (111) 1.0 0.8 0.6 Interface free energy [J m -2 ] 0.4 0.2 0.0 -1.5 -1.0 -0.5 ( 3 3 4 ) (4 34 ) (

(

(

(414 liquid structure to the same level of importance as the crystal surface morphology in understanding crystallization, a knowledge that can enable material advances through the incorporation of liquid-structure engineering as a novel pathway for synthesis. Our results were only made possible by employing atomistic simulations and ML together. The strength of this combined approach is that one can perform complex simulations and yet glean physical insight from notoriously haphazard atomic environments. This innovative application of ML in materials science blends conventional scientific methods with data science tools to produce physically consistent predictive models and novel conceptual knowledge.

Methods

Silicon crystal growth simulations. The MD simulations were performed using the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS 35 ) software, with the interactions between silicon atoms described by the Stillinger-Weber 36 interatomic potential. The timestep was selected as ~1/56th of the period of the highest-frequency phonon mode of this system, or Δt = 1 fs. The crystal growth simulations contained 500,000 atoms and were initialized with a spherical crystalline seed of ~3000 atoms in the diamond cubic structure. The lattice parameter for the atoms in the crystal seed was chosen taking into account dilation due to thermal expansion, then the remainder of the simulation cell was filled with randomly distributed atoms at the equilibrium liquid density for that temperature at zero pressure.

The system was equilibrated by first relaxing the liquid atoms using a Conjugate Gradient 37 algorithm for 200 steps. Next, the liquid was equilibrated at finite temperature using the Bussi-Donadio-Parrinello 38 (BDP) thermostat for 3 ps with a damping parameter of 0.1 ps. Finally, liquid atoms were equilibrated for 2 ps at zero pressure and finite temperature using the same thermostat just described and a chain Nosé-Hoover barostat [39][40][41][42][43] with damping parameter of 1 ps and a chain length of three, allowing only for isotropic dilation/contraction of the system. During the entirety of this equilibration process the crystalline seed atoms were kept frozen at their equilibrium crystal structure with fixed lattice parameters (i.e. they did not dilate/contract with the liquid atoms). After equilibration the BDP thermostat and chain Nosé-Hoover barostat were applied to the entire system, both with damping parameter of 1.0 ps, to maintain the system at finite temperature and zero pressure for a total of 3 ns during which snapshots were recorded every 1 ps. The crystal growth process can be seen in the Supplementary Video 1. Snapshots saved from the MD simulations were subsequently relaxed using 20 steps of the Steepest-Descent 37 algorithm. The crystal growth simulations were performed at temperatures ranging from 1125 to 1500 K in intervals of 25 K.

The damping parameter for the thermostat was selected conservatively such that the liquid diffusivity was not affected by the presence of the thermostat, i.e. it had the same value within the statistical uncertainty as the diffusivity computed without a thermostat. Thus, the thermostating of the crystal growth simulation was performed gently as to not affect the kinetics of the system. See Supplementary Note 6 for more details on the diffusivity calculations.

Phase identification. In order to identify to which phase (liquid or crystal) each particle belongs, we used the order parameter introduced by Rein ten Wolde et al. 44 . The complete description and analysis of the construction of this order parameter can be found in the Supplementary Note 2. Ultimately, this method provides us with a parameter α i (t) for each atom i at the time t of each MD snapshot. The physical interpretation of this parameter is that α i is the fraction of bonds that atom i makes that resemble bonds in a perfect crystal structure. As shown in the Supplementary Note 2, the parameter α i correctly identifies atoms in the perfect crystal or bulk liquid with accuracy of 100% within the statistical uncertainty. It is important to notice that although α i can discern between liquid and crystal atoms, it does not differentiate between crystalline structures. We confirm that the silicon atoms are indeed crystallizing in the diamond cubic structure by performing the polyhedral templatematching 45,46 analysis.

Encoding atomic dynamics (ML labeling). The dynamics of each atom was encoded using the time evolution of the α i (t) order parameter. A representative plot 47 of α i (t) is shown in Supplementary Fig. 4b and c. Notice that due to thermal fluctuations the instantaneous value of α i (t) for atoms in the liquid and crystal phases might differ from their perfect values of 0.0 and 1.0, respectively, even after the short Steepest-Descent relaxation. Hence, we perform a moving-window average of α i (t) with window length of 20 ps and use the window-averaged α i ðtÞ to label the atomic dynamics as illustrated in Supplementary Fig. 4a. Atoms with α i ðtÞ ¼ 0 for t ∈ [t 0 -τ ℓ , t 0 + τ ℓ ] receive label y i = -1 at time t 0 . These are atoms deep in the liquid phase that will not be transitioning to the crystal state in the near future, neither have tried to transition in the near past. From the analysis of curves such as in Supplementary Fig. 4c we choose τ ℓ = 15 ps as a reasonable value. Next we identify atoms that have just started to move out of the bulk liquid (i.e. crystallizing atoms) as those within a 20 ps window from the point where α i ¼ 0:25, i.e. y i = 1 for t ∈ [t 0 -τ, t 0 + τ] where α i ðt 0 Þ ¼ 0:25 and τ = 10 ps. See Supplementary Note 2 for more details on the labeling process.

Local-structure fingerprint (ML features). The local structure surrounding each atom was characterized using a set of radial structure functions 13,16 :

where i is the atom whose local structure is being described, n(i) is the number neighbors of i within a cutoff radius r cut , r ij is the distance between atom i and one of its neighbors j, r and σ are two parameters that define the radial structure function. These smoothly varying functions of r count the number of neighbors of i at a distance r. In this interpretation, parameter r represents the radial distance from i at which we are counting the number of neighbors, while σ adjusts how smoothly the function varies as atoms move in and out of the distance r vicinity.

We have used a grid-search 48 algorithm to perform the hyperparameters tuning (see Supplementary Note 1 for more details), resulting in σ = 0.5 Å, r cut = 10.8 Å, and r n = (2.0 + 0.4n) Å, with n = 0, 1, 2, …, 20. With this set of 21 radial structure functions-one for each value of r-the local-structure fingerprint of each atom i was built as a vector:

Softness calculation. The data was assembled by pairing the dynamic labels y i with their corresponding structural fingerprint x i . Then, 10,000 ðy i ; x i Þ pairs were randomly selected and equally divided between the y = -1 and y = 1 classes to train a support vector machine [17][18][19] (SVM) classifier. The SVM algorithm finds the hyperplane of the form w Á x À b ¼ 0 that optimally separates the two classes, where w and b are the parameters that define this hyperplane. Before training the SVM classifier the elements of the fingerprints were standardized 48 to have zero mean and standard deviation of one. The optimal hyperplane found, denoted by the parameters w * and b * , correctly separates the two classes with an accuracy of 96%. See Supplementary Note 1 for more details about how these parameters are found and an in-depth analysis of the quality of the classifier computed. All results shown here were obtained using data from the crystal growth simulation at T = 1500 K to train the SVM classifier. However, the results can be reproduced within the statistical uncertainty when training at any other temperature, as shown in the Supplementary Note 3.

Once the SVM classifier has been trained it was applied to the entire data set, composed of 27.5 million data points (excluding the data used for training, hyperparameter tuning, and cross validation). The value of softness 14 for each data point (or atom) is the signed distance from the hyperplane, or S i ¼ w à Á x i À b à for each atom i.

Parameter estimation and temperature extrapolation. In order to test how predictive the LSD and WF models are we performed the model parameterization of both models using only the data collected for low undercooling (i.e. T ≥ 1388 K) and observed how well the model predicts the temperature dependence for higher undercoolings (i.e. temperatures as low as 1128 K). As shown in Fig. 3a the LSD model is capable of predicting the growth rate at temperatures it was not parametrized on, while the WF model only reproduces the simulation results at temperatures it was fitted to. We attribute this to the fact that the LSD model accounts for the Arrhenius family of thermally activated atomic events leading to crystal growth, as labeled by S. This is fundamental physical information that is not incorporated in the WF model.

For Figs. 5 and 6a only, the S dependence of the parameters of the LSD model (i.e. ΔE a , k 0 , and γ) was measured after reparametrizing this model using the data from all simulations with T ≥ 1206 K. The reparametrization was necessary only in order to reduce the statistical uncertainty of the parameters measured. This range of temperatures was chosen because below T = 1206 K the kinetic factor kðT; SÞ showed signs of non-Arrhenius behavior for some values of S due to the approaching glass transition temperature. Notice that below T = 1206 K the bulk liquid diffusivity also shows signs of departure from the Arrhenius behavior. Thus, there are no reasons to expect the Arrhenius behavior for kðT; SÞ to hold below T = 1206 K because the mobility of liquid atoms close to the crystal surface is smaller than in the bulk due to IIO effects.

Principal component analysis. Figure 4b was obtained by applying the PCA 48 to a data set containing equal amounts of crystallizing, liquid, and crystal, for a total of 60,00 data points. Crystal atoms were defined as those for which α i ðtÞ ¼ 1:0 (bulk crystal atoms) or α i ðtÞ ¼ 0:75 (stacking fault atoms) for t ∈ [t 0 -τ x , t 0 + τ x ] with τ x = 15 ps. From the PCA we obtained the components of each data point along the two eigenvectors with largest eigenvalues, which are used to plot Fig. 4b. The atom trajectory was obtained by applying the same PCA transformation along a single 3 ns trajectory of an atom in a simulation at 1500 K.

Growth rate determination. The number of atoms in the crystallite N(t) at any given time t was determined as the number of atoms with α i ðtÞ>0:25. From this information the effective crystallite radius R eff ðtÞ (shown in Supplementary Fig. 10a) was computed assuming a spherical shape (which results in κ ¼ 2=R eff , where κ is the geometrical factor in the Gibbs-Thomson term). The growth rate for each temperature, shown in Fig. 3a, was determined by a linear fit of N(t) over the time interval for which R eff 2 ½80A; 100A. This interval is such that the crystallite is small enough to not be affected by finite-size effects, but large enough to give the system time to equilibrate into a steady-state growth condition. See Supplementary Note 4 for a more detailed analysis.

Error bars in Fig. 3a represent the 95% confidence interval as computed using the bootstrap method with 1000 samples of the same size as the original distribution.

Interface temperature. When studying crystal growth, it is important to differentiate between the temperature of the supercooled liquid surrounding the crystal (but far from the interface) from the solid-liquid interface temperature. Under steady-growth conditions these two temperatures will differ because of the latent heat released at the interface and the finite rate of heat transport. Here, the interface temperature was computed by considering only the kinetic energy of atoms with α i 2 ð0:15; 0:75Þ. This interval of α was selected because it includes both, interfacial liquid and interfacial crystal atoms. Supplementary Fig. 10b shows the interface temperature under steady-state growth as a function of the surrounding liquid bath temperature. See Supplementary Note 4 for more details. The crystallizing atom (green) has its local structure (illustrated here only by its first neighbors) affected by the nearby solid-liquid interface. This effect is anisotropic: high-index or rough interfaces interact strongly with the liquid and cause significant ordering of the liquid, which becomes rigid, resulting in large activation energies ΔE a when compared to the barrier for diffusion in the liquid ΔE d . b Low-index interfaces interact weakly with the liquid and cause very small ordering, resulting in low ΔE a . c Schematic illustration of the effect of interface-induced ordering of the liquid on the free-energy landscape of crystallization. Note how the liquid free-energy basin moves to the left with increasing ordering, causing ΔE a to become progressively and continuously larger.

Crystal surface analysis.

3.0 Å and a smoothing level of 10. From this mesh the surface directions were inferred and averaged over the time interval for which the crystal growth occurs in a steady state. The data for constructing Fig. 6b was obtained by finding the orientation of the closest surface to each crystallizing atom. The only atomically smooth surfaces in silicon are {111} surfaces 22 , consequently steps can only exist in these surfaces. Hence, the step-step separation distance shown in Fig. 6c was computed assuming that (111) faceting occurs at all surface orientations. Solid and liquid free energies. The accurate calculation of the solid and liquid free energies is crucial in crystal growth studies. As shown in the Supplementary Note 5, employing approximations such as the quasi-harmonic approximation results in the underestimation of the predicted growth rates by as much as 36%. For this reason, we performed the solid and liquid free energies calculations using state-ofthe-art nonequilibrium thermodynamic integration methods that make no approximating assumptions on the physical characteristics of the system. The crystal free energy was determined using the nonequilibrium Frenkel-Ladd [50][51][52] (FL) and the reversible scaling 52,53 (RS) methods, following closely the approach described by Freitas 52 . For both methods a system of 21,952 silicon atoms in the diamond cubic structure was employed. The thermodynamic switching was performed in 200 ps for each direction, before which the system was equilibrated for 20 ps. The FL switching was realized for temperatures ranging from 100 to 2000 K in intervals of 100 K. For each temperature the switching was repeated in five independent simulations to estimate the statistical uncertainty. Similarly, the RS switching was also repeated five times. The S-shaped function was employed in the FL switching, while the RS switching was performed with T i = 100 K and T f = 2000 K under the constant dT/dt constraint. The system's center-of-mass was kept fixed for the FL and RS simulations, while a Langevin 54,55 thermostat with damping parameter of 0.1 ps was applied. For the RS method a chain Nosé-Hoover barostat with damping parameter of 1 ps and chain length of three was used to keep zero pressure. The absolute free energies and a comparison with the harmonic and quasi-harmonic approximations 56 can be seen in Supplementary Fig. 11.

Liquid free energies were computed using the Uhlenbeck-Ford 57,58 (UF) and RS 53,58 methods, following closely the approach described by Leite 58 . The liquid free energy calculations had the same number of atoms, switching time, equilibration time, and thermostat as the crystal free-energy calculations. The liquid density was the equilibrium density at zero pressure, with the thermal expansion dilation taken into account. For the UF method we used p = 50, σ = 1.5 Å, and a cutoff radius of r c = 5σ. The UF switching was performed linearly with time, while the RS switching had the same time dependence as the crystal with T i = 2000 K and T f = 1100 K (the lower final temperature T f was chosen to avoid the liquid vitrification at low temperatures). For both methods-UF and RS-the switchings were repeated in five independent simulations to estimate the statistical uncertainty.

Copper. All results for copper were obtained from simulations that followed the exact same specifications as the simulations described above for silicon. The only modifications performed are described in this section.

The interaction between copper atoms was described using the embedded-atom method 59 interatomic potential of Foiles et al. 60 . The timestep was selected as approximately 1/66th of the period of the highest-frequency phonon mode of this system, or Δt = 2 fs. The crystal growth simulations contained 1,000,000 atoms (notice that this is twice the size of the simulations for silicon) and were initialized with a spherical crystalline seed of ~24,000 atoms (eight times larger than for silicon) in the face-centered cubic structure. The difference in system size allowed us to explore much lower undercoolings for the crystal growth simulations, which ran for a total of 2 ns per temperature. In contrast to the simulations for silicon, the snapshots saved from MD simulations were not subsequently relaxed before computing structural parameters because it has been shown that energy minimizations lead to significant crystallization in metallic systems 61 . The crystal growth simulations were carried out at temperatures ranging from 900 to 1200 K in intervals of 25 K. The growth rate for each temperature, shown in Fig. 7a, was determined by a linear fit of N(t) over the time interval for which R eff 2 ½100Å; 120Å, as detailed in Supplementary Note 4.

The local-structure fingerprint x i for the copper atoms was composed of a set of 37 radial structure functions defined by the following parameters: σ = 0.3 Å, r cut = 20.6 Å, and r n = (2.0 + 0.5n) Å, with n = 0, 1, 2, …, 37. These parameters were obtained through the same hyperparameter optimization process applied to silicon, as described in the Supplementary Note 1. The SVM classifier trained with the data collected for copper had accuracy of 97%. The results presented here were obtained using data from the simulation at T = 1200 K to train the SVM classifier. The total size of the data set for copper was 77.9 million data points.

All results for the LSD model for copper (including Fig. 7 and Supplementary Fig. 19c) were obtained using data from simulations at T ≥ 1136 K. Below this temperature the kinetic factor showed signs of non-Arrhenius behavior for some values of S due to the approaching glass transition temperature, similarly to what was observed for silicon. Notice that below T = 1136 K the bulk liquid diffusivity also shows signs of departure from the Arrhenius behavior. Thus, there are no reasons to expect the Arrhenius behavior for kðT; SÞ to hold below T = 1136 K because the mobility of liquid atoms close to the crystal surface is smaller than in the bulk due to IIO effects.

The solid and liquid free energies were computed for systems containing 19,652 copper atoms. The FL and RS methods were applied with a switching time of 400 ps for each direction, preceded by an equilibration time of 40 ps. The FL switching was realized for temperatures ranging from 100 to 1300 K in intervals of 100 K. The RS switching for the solid was performed with T i = 100 K and T f = 1300 K, while for the liquid we used T i = 2000 K and T f = 900 K. For the UF free-energy calculations we used p = 75, σ = 1.3 Å, and a cutoff radius of r c = 5σ.