Nucleation pathways in barium silicate glasses

Nucleation is generally viewed as a structural fluctuation that passes a critical size to eventually become a stable emerging new phase. However, this concept leaves out many details, such as changes in cluster composition and competing pathways to the new phase. In this work, both experimental and computer modeling studies are used to understand the cluster composition and pathways. Monte Carlo and molecular dynamics approaches are used to analyze the thermodynamic and kinetic contributions to the nucleation landscape in barium silicate glasses. Experimental techniques examine the resulting polycrystals that form. Both the modeling and experimental data indicate that a silica rich core plays a dominant role in the nucleation process.

Nucleation and crystallization in glasses is a rich subject area covering tunable materials for high technology applications to fundamental scientific inquires [1][2][3][4][5][6] . The final nucleated glass, or glass-ceramic, may have improved material properties, such as toughness 6 , over the parent glass. By understanding the nucleation event, it may be possible to design an improved glass-ceramic or recognize how the devitrification process occurs, which is crucial for the glass making industry. The Classical Nucleation Theory (CNT) 1,7-9 is most commonly used to describe the nucleation and growth of the emerging new phases. Qualitatively, CNT provides a good framework to understand the process. However, problems arise when this method is used to obtain a quantitative description (thermodynamic and kinetic) of the nucleation phenomenon. In this work, we focus on homogeneous nucleation rates that occur in the glass far below the melting temperature (∆T = T m − T, where T m is the melting temperature). At these temperatures, the critical nuclei are only 1-2 nm in size, making it difficult to experimentally probe their size evolution.
According to CNT, the thermodynamic properties of the evolving new crystalline phase are considered size independent and equivalent to the properties of the macroscopic phase, with a sharp interface separating the nucleating cluster from the parent phase (i.e., the capillarity approximation). For instance, in CNT the interfacial free energy between a crystal cluster and the glass will be the same for very small clusters as it is for macroscopic ones. It is typical to express the properties of a cluster in terms of its radius, r. The work of cluster formation W(r), for a spherical cluster and a constant interfacial free energy term σ , as shown in Eq. (1).
The driving free energy for nucleation and growth can be computed using the difference of the Gibbs free energies of the bulk crystal and the amorphous phase, g . The critical cluster radius, r * , are cluster sizes below which they tend to dissolve and above which they are grow.
The assumption of a sharp interface between the cluster and the parent phase is a reasonable one for the case of vapor-liquid nucleation that CNT was originally developed to address. However, in dense phase systems where the parent phase is in intimate contact with the developing new phase, the interface cannot be as sharp. This was originally proposed by Turnbull to explain differences between his experimental results and predictions of the CNT for nucleation in liquid mercury 10 . To better describe the emerging cluster and interface region, two theories have been developed: the Diffuse Interface Theory (DIT) and Density Functional Theories (DFT). DIT does not assume the capillarity approximation, but instead describes the work in terms of an assumed profile of the Gibbs free energy as a function of distance from the cluster center [11][12][13] . The DFT theories [14][15][16] are based on an order parameter description of the phase transition. Both DIT and DFT predict a diffuse interface. We chose DIT since it could be more readily be computed.
A long-standing problem of using CNT for glass-ceramic systems is the work of critical cluster formation, W(r), as a function of temperature. Based on experimental studies at low temperatures (lower than the maximum www.nature.com/scientificreports/ nucleation temperature, T max ), W(r) begins to increase with decreasing temperature, rather than continuing to decrease as expected [17][18][19][20][21] . This discrepancy has had at least five different possible explanations. These include: (1) a change in the sign of the temperature dependence for the interfacial free energy with a temperature dependence 10 , (2) a change in the driving free energy with temperature 20 , (3) a cross-over in the kinetics of long-range diffusion and interfacial attachment 22 , (4) a decrease in the thermodynamic driving free energy due to elastic stress 20 , and (5) a variation of the cooperatively rearranging regions of the critical nuclei that may account for kinetics and longer annealing times 23 . However, these explanations alone have not been able to clarify the description of nucleation at such deep undercoolings (50 degrees below T max ). A common thread in these explanations is the thermodynamic landscape as a function of temperature, which is explored in this work. By illuminating the thermodynamic nucleation contributions as a function of chemistry and temperature, one may be able to detect differences in the nucleating pathways for the sub-critical nuclei. Considering a pathway change, the growing cluster composition is important as it may elucidate other routes (i.e., micro-phase separation leading to nucleation). Over recent years, a "two-step" nucleation mechanism has been gaining attention, where a disordered droplet will form first, followed by crystallization 4,[24][25][26] . Probing these intermediate clusters is important in identifying these paths. How distinct are these two steps of nucleation? Can they occur concurrently? Most nucleation theories assume that the clusters grow by monomer attachment and detachment. This is a straightforward counting method for the nucleation of droplets of water in vapor, where the monomer is well defined. But, how well does this scheme work in glass-crystal nucleation? Does one count neutrally charged oxide units (every i-SiO 2 or i-BaO); does one track individual elements; or would one count by unit-cells (every i-{5BaO 8SiO 2 } or i-{1BaO 2SiO 2 } units)? The barium silicate system is known to contain multiple crystalline polymorphs within the same crystal 27 , so how would monomer tracking look in this case? How would the evolution of the nuclei modify the interfacial region and resulting instantaneous interfacial free energy? To best explore this nucleation cluster free energy landscape, we choose to track the cluster evolution by every oxide unit. Between CNT and DIT, experiment and models, we will probe the nucleation pathways. Reality is more complicated and multiple concurrent intermediate phases will be shown to play a role when the final state is reached 4,24 . Classical nucleation theory background. Gibbs 28,29 studied phase transitions that occur under near equilibrium conditions, where the probability of fluctuations that lead to the new stable phase is small. This corresponds to a large nucleation barrier. He found that the barrier decreases when the system is quenched more deeply into the metastable region. Later, Volmer and Weber 30 constructed CNT based on Gibb's reversible work for a critical cluster and a new kinetic model that they introduced. In this formulation, the work of formation, W(n), of a cluster of the new phase containing n particles, is used to describe an equilibrium cluster size distribution for clusters smaller than the critical size. The critical cluster size, n*, or denoted by the cluster radius, r*, corresponds to the maximum work of formation W* (the asterisk in the superscript denote a critical size value). Assuming a steady-state cluster distribution, Becker and Döring 31 showed that the steady-state nucleation rate, I, is proportional to the Boltzmann weighted probability of the critical fluctuation and a kinetic prefactor, A*, where T is temperature and k B is the Boltzmann constant. Turnbull and Fisher extended this nucleation formulation to include crystallization from a liquid by relating the prefactor to kinetic processes in the liquid 32 . For deep undercooling, the influence of the kinetic prefactor is difficult to assess. The key kinetic term in the prefactor is the diffusion rate, which for crystallization from liquids and glasses is typically determined from the viscosity, assuming the Stokes-Einstein relation. However, the viscosity is difficult to measure and model at low temperature. Also, it is believed that there is a breakdown of the Stokes-Einstein relation at these temperatures 33,34 . Diffuse interface theory background. Based on experimental nucleation studies, Turnbull argued that instead of having a sharp interface between the parent and new phase, as assumed by Gibbs, the interface was actually somewhat diffuse. Gránásy 11,12 and Spaepen 13 independently proposed a phenomenological model to account for the changing thermodynamic quantities through the diffuse interface. Here, the work of cluster formation, W, is a function of the Gibbs free energy, g, The Gibbs free energy at radius r is g(r) = ∆h − T∆s, where ∆h and ∆s are the enthalpy and the entropy difference between the glass and crystal respectively. Within the DIT, there is an offset in the thermodynamic quantities, ∆h and ∆s, as a function of distance from the crystal cluster center, r. While the precise profile for g(r) through the interface is unknown, by approximating ∆h(r) and ∆s(r) as a series of step functions Eq. (3) can be directly solved, assuming that the Gibbs free energy between the bulk glass and the center of the crystal is known (i.e. from the Turnbull approximation or from specific heat data of the glass and crystal).

Results
This section will discuss (1) the results of the GCMC model for barium silicate nucleation, (2) a comparison to the DIT and ordering between the solid and liquid phases, (3) low temperature concerns, and (4) the experimental results. Different polymorphs of the barium silicate glass are considered: 1BaO·2SiO 2 , 3BaO·5SiO 2 , 5BaO·8SiO 2 , and 4BaO·6SiO 2 . For convenience, these different crystals will be referenced by their barium silicate stoichiometry. For example, 5BaO·8SiO 2 will be labeled 5-8. www.nature.com/scientificreports/ Barium silicate GCMC study. Nucleation modeling is sensitive to the force field used. It was found that the melting temperature for the 1-2 was 60 K higher than the experimental value (1693 K). This is reflected in the nucleation temperatures in the model, where it is predicted that the maximum nucleation temperature, T max , occurs at 1090 K, higher than the experimental value (985 K 21 ). Following CNT it was found that the maximum work per crystal polymorph occurred between 2 to 3 unit cells (4 to 6 unit cells for the smallest 1-2 system). The maximum work of formation (i.e. the work at the critical size) for the 1-2 and 5-8 glasses in shown in Fig. 1. This is in good agreement with previous experimental results 21 , where the maximum work at 1090 K for the 1-2 and 5-8 glasses are 32 k B T (vs. the experimental value of 31 k B T) and 24 k B T (28 k B T experiment) respectively. Interestingly, from Fig. 1, there is a continued decrease in the work of formation. This lower temperature had to be sampled 30% longer than the other temperatures due to slow convergence. This was caused by a slowly evolving cluster structure with very cohesive elements (i.e., the force field is too attractive). While the overall melt stoichiometry may be that of a barium silicate crystal, it is heterogeneous and there is a possibility that a region may contain a higher than average concentration of barium (or silica). As shown in Fig. 2, this will guide where and how the nucleation occurs. Note that W* in Fig. 1 is the work of the cluster (Eqs. 7-9) at the first maximum in the Nucleation Free Energy (NFE) contour map shown in Fig. 2. The NFE map is created by thermodynamic integration of the cluster's work. This Fig. 2 only has the thermodynamic contributions as there are no continuum diffusion kinetics in the GCMC model. According to CNT, one must follow the diagonal line representing the polymorph stoichiometry. For all polymorphs and at all three temperatures, there is a protected valley labelled as A. The first unit cell of each of these polymorphs occur in this region (except, it is 4 unit cells of the 1-2 system). These pre-critical clusters then must climb the NFE barrier region (marked with B in Fig. 2) to reach more stable cluster configurations.
Comparing the 980, 1090, and 1160 K contour plots in Fig. 2, the NFE landscape becomes smoother with increasing temperature, and the size of the higher energy regions (red contours) decreases. A decreasing NFE barrier at higher temperatures will cause the nucleation rate to increase. However, as the CNT polymorph stoichiometry lines are traversed across Fig. 2, different paths may be taken to reach the next cluster size (e.g., first add SiO 2 then BaO or vice versa). No matter which route is taken to go beyond this large NFE barrier in the 980 K case, it is necessary to either travel through it or dissolve BaO units to form a more silica rich cluster. With increasing temperature, more routes become accessible since the NFE contours decrease.
Beyond this barrier at B, the clusters will reach another NFE minimum at C. This blue region extends diagonally towards the silica axis, marked by D. The system should follow the lowest free energy path. With such a low NFE region at D, a glass heterogeneity of this composition, which is quite distinct from that of the bulk, may in fact serve as the first step in the nucleation process. Given a specific bulk barium silicate stoichiometry, is it possible to have regions that are enriched in silica and deficient in barium? Glass melting research has been working on reducing compositional inhomogeneities for a long time and this remains a concern for glass manufacturing [35][36][37][38][39][40][41] . However, it is possible to find such regions within the glass structure. Depending on the raw material, how the glass is melted (re-melted), and the annealing cycles used, all may affect the compositional inhomogeneity. These factors would also affect how the cooperatively rearranging regions in the glass melt evolve 42,43 . Analyzing compositional variation in glass melts, however, is outside the scope of this work, though is a critical aspect of our ongoing experimental studies of nucleation in the barium silicate system.
The ordered structure should have comparable amounts of baria (BaO) and silica (SiO 2 ) if it grows according to CNT throughout the cluster volume. Here, we analyze the last 1000 Monte Carlo configurations by tracking the amount of silica to the total sum of silica and baria found in the inner 50% of the cluster volume. The standard deviation of this analysis is less than the resulting contour intervals. For clusters smaller than 10 silica and 10 baria the inner volume is very small and numerical artifacts are seen. Surprisingly, for all three www.nature.com/scientificreports/ temperatures shown, the inner volume is silica rich (Fig. 3). Using the diagonal line as the 5-8 unit cell line (the other polymorphs would fall nearby on either side), it appears nucleation seems to be mainly driven by a silica rich environment. This compares favorably with the study by Avramov et al. 44 on topological constraint theory 45,46 and nucleation where they found that the network needs to be slightly rigid to enable crystallization. The silica structures would be more rigid than the barium oxide network modifier structures, simply due to a higher degree of network connectivity provided by fewer non-bridging oxygens and more fully polymerized silicate tetrahedra. Additionally, the formation of silica (e.g. cristobalite) would be the most favored structure on this NFE landscape ( Fig. 2 labelled D). Finding cristobalite or quartz during the glass making process is common and undesired (a.k.a. devitrification).
Comparison to DIT and ordering between phases. Clearly, CNT has its deficiencies, especially in describing the atomic environment of the pre-critical clusters. DIT brings in a more detailed thermodynamic description of the order parameter as a function of distance through the cluster and into the melt. Here, we used a molecular dynamics simulation where we explicitly model both the crystal cluster and the melt (Fig. 4). In DIT the critical cluster size is determined in the usual way, by setting the derivative of the work of cluster formation at a given temperature equal to zero and solving for the radius (or number of particles/monomers). Here, DIT predicts that the critical cluster size has a radius ≈ 12 Å at 980 K, which is slightly larger than the biggest clusters modeled with the GCMC approach (radius ~ 10 Å). The critical size is smaller from the GCMC approach compared to DIT. It is reasonable considering the NFE landscape has hills and valleys (fluctuations) that will lead to a critical nucleus. The same fluctuations were seen in a previous lithium silicate study 26 as the cluster grows. Additionally, the NFE landscape in Fig. 2 seems to become smoother at larger sizes. This qualitatively compares well to the MD potential energy results in Fig. 5 (discussed below). To further illuminate these fluctuations and the developing crystal-glass interface, we studied a critical size nuclei (r = 12 Å) and a post-critical one (r = 41 Å).
Using the MD glass-ceramic model, we radially averaged the potential energy and bond-orientational order parameter Q4 (Fig. 5). Since this is an average over a spherical radius, depending on the bin it may cross over different crystal planes-the Q4 parameter is especially sensitive to this. The Potential Energy of a crystal is equivalent to the enthalpy, assuming that no pressure-volume work is being performed. For DIT's critical cluster of 12 Å, there are highs and lows of energy corresponding with crystal lattice positions. The magnitude and pattern clarity continue into the melt (right side of dotted line in Fig. 5a). There is a decrease in the Q4 structure when passing through this interface (Fig. 5b). The nearest melt layer shows signs of the Q4 and Potential Energy patterning of the crystal. This shows that within the simulation, the cluster is stable and the MD model agrees with DIT's prediction. Interestingly, there is little difference between the crystal and the melt in Fig. 5a,b (left vs. right sides), meaning that the cluster is not too dissimilar to the glass. This indicates that the cluster is more like the interface than a crystalline cluster with a sharp interface. It is worth mentioning these values extend into the melt, inferring a larger interface, or fluctuation.
The post critical cluster (radius = 41 Å, Fig. 5c,d) shows oscillations of the Potential Energy that follow the crystalline ordering; the magnitude is lower in the melt (being randomly distributed). Additionally, the larger crystal cluster shows a clear loss of order. This ordering begins to fall off around 30 Å. In other words, the cluster has an interfacial width of approximately 10 Å! As the cluster grows, the interface becomes sharper and less significant relative to the cluster radius, being more consistent with the assumptions of CNT. This finding is in line with previous studies 1,2 . Interestingly, the crystal's oscillating Potential Energy extends into the melt (~ 10 Å, Fig. 5c) further than the bond order does. This may indicate that the correct layering occurs first, followed by an increasing density (lowering of entropy). Then, the crystalline ordering occurs to reduce the total volume. This is opposite of Tanaka's 3 work, where the crystallization and growth occur with a density fluctuation followed by the crystalline bond ordering. However, this work agrees with a metallic liquid nucleation study 4 where an icosahedral ordering occurs first, followed by crystallization. Clearly, there are nucleation nuances that remain open.
Lower temperature concerns. From Fig. 5, there is no temperature effect on the width of interfacial region in the atomistic simulation. One might expect the interfacial width to widen as temperature increases. Additionally, the lowest temperature in the GCMC study had slow convergence since it samples a flatter configuration energy space. To better address this concern, a simple and straightforward investigation of the barium silicate diffusion coefficient allows one to assess the force field at lower temperatures and gives insight into the nucleation kinetics.
The Ba diffusion is shown in Fig. 6; Ba diffuses faster than O and Si. To obtain the diffusion coefficient a slope of unity must be reached on a plot of Mean Squared Displacement (MSD) vs time. If the slope is less than unity the atoms are trapped in a cage and the diffusion coefficient cannot be measured until they escape. At only the two highest temperatures has the slope been attained for a proper calculation of the diffusion coefficient. The lines in Fig. 6 are composed of three restarted simulations (of 200 ns each). If one approximates a linear relationship as the slope changes from each additional run, for the T = 1190 K case it would require an additional 4600 ns trajectory to reach a slope of unity. That is assuming a linear relationship, which is unlikely because these cooler temperatures become nearer to the glass transition temperature, it is more of a measure of the bare minimum time needed. It was found that the diffusion coefficient for Ba was 0.01378 nm 2 /ns and 0.14563 nm 2 /ns at temperature of 1500 and 1773 K, respectively. This value is two to four orders of magnitude faster than an estimation made using the crystal growth velocity 5 , where it was extrapolated to be 7.94 × 10 -4 nm 2 /ns at T = 1693 K. Clearly, the modeled diffusion coefficient is faster than the experiment suggests. Predicting lower temperature diffusion would suffer from similar inaccuracies and differences in the glass transition temperature. The difficulties in the GCMC thermodynamic model can be effectively mitigated with some additional sampling, however the dynamic contributions would be orders of magnitude more time consuming.     www.nature.com/scientificreports/ Experimental study: 5-8 system. Shown in Fig. 7 are SEM images of the 5-8 composition after a nucleation treatment for 2 h at 780 °C. These images compare the morphology with and without acid etching, showing large changes when briefly exposed to the 0.1% Hydrofluoric acid (HF) solution (etchant). There appears to be a center region in the fractured and not etched sample (Fig. 7). The acid etching generates voids and sharpens the resolution of the nuclei center region versus the outer region of the spheroidal regions. Silica is especially susceptible to HF acid, suggesting that an explanation for these differences is that the cores of these nuclei are Si-rich and thus more readily removed during etching. To look for such compositional differences, additional SEM and EDS analyses were conducted on the nucleated 5-8 glass, the results of which are given in Fig. 8. Indeed, the EDS analysis of the center region versus the outer parts of the spheroids confirms a significant compositional difference, with the core region being enriched in Si and depleted in Ba. This is certainly consistent with the observation of different etching behavior across the nucleus. Further, the results of the GCMC study also found that all clusters at any temperature contain a silica rich core (Fig. 3) and that a silica rich region might serve as a better nucleation site (Fig. 2, label D). There is a size discrepancy between the model and the SEM image, where the modeled cluster radii towards the upper left corner in Fig. 2 are ~ 12-20 Å, while the radii of the clusters shown in Figs. 7 and 8 are 600 nm. In spite of these size differences, good agreement between the predictions of the model and the experimental data is observed for compositional fluctuations within the 5-8 nuclei. www.nature.com/scientificreports/ Experimental study: 1-2 system. The 1-2 system was nucleated at 850 °C for 30 min. Other samples were also nucleated at lower temperatures and longer times and showed similar crystalline features. The morphology of the 1-2 crystal is needle-like and distinct from that of the 5-8's spheroid shape. These needle-like crystals also exhibit the same core and shell organization as the 5-8 (Fig. 9a). A closer inspection using Bright Field STEM shows a variety of crystal types and twinning planes. The 1-2 unit cell line in Fig. 2 is further away from the other barium silicate crystal lines in the same figure. Yet, a similar polycrystalline sample emerges ( Fig. 9) in that we find needle-like and spherulite crystals having 1-2, 5-8, and 2-3 compositions. This occurs due to evolving local composition (and chemical potential) and the ability to transverse to nearby CNT crystal polymorph pathways (i.e., there are no large barriers that are parallel to the different barium silicate crystal line in Fig. 3).

Discussion and future work
In this paper, we demonstrate the richness of nucleation by: (1) showing how a nuclei develops through CNT, DIT, different modeling efforts, and experiment; (2) identifying the feasibility of an alternative silica rich pathway; (3) identifying different nearby CNT pathways leading to co-crystallization of other Ba-Si forms, and (4) demonstrating that the thermodynamic contributions to nucleation also decrease and cannot describe the apparent T max anomaly. The Classical Nucleation Theory discretizes the macroscopic world into a microscopic realm, but, these results show that nucleation is a complex process having a greater depth than CNT ever envisioned. The discretization is, in fact, a spectrum of possible microstates all interacting with one another, as shown in the GCMC model. The Diffuse Interface Theory begins to fill in the spectrum by adding an interfacial region and connecting the emerging cluster to the melt; as the cluster grows the interfacial width decreases. The modeling shows that for sub-critical and critical cluster sizes, the interface is a significant portion of the cluster. The modeled work of cluster formation as a function of cluster composition can lead to a variety of nucleation pathways, where each cluster has a different thermodynamic driving force and interfacial free energy, all dependent on the local composition. This is shown from the experimental investigation, where independently from overall melt stoichiometry, other barium silicate crystals are identified. An interesting behavior, common between the different pathways, is the silica rich core. Experimental measurements have shown that the inner core of the 5-8 cluster is silica rich, which qualitatively agrees with predictions from the GCMC model. If the CNT discretized lines are not followed, but the whole barium silicate landscape is followed, the low NFE region on the silica rich axis, having low barium content (Fig. 2, label D), may serve as the first step in the nucleation process. Qualitatively, this agrees well with studies of nucleation in supercooled water, where immobile regions allow water molecules to nucleate. In this case, the silica rich regions of a melt are more thermodynamically stable and less mobile (being a network former), allowing barium silicate to nucleate from them. A scientific question then arises, being a silica rich heterogeneity, would this be a result of a poorly mixed glass occurring via heterogeneous nucleation instead? This may be the guiding principle of cluster growth and devitrification.
The number of different pathways as a function of temperature and composition tends to decrease as a function of temperature, shown by the more pitted and hillier NFE landscape in Fig. 2 as temperature decreases. www.nature.com/scientificreports/ However, the work of formation for the critical cluster does tend to decrease with decreasing temperature, unlike what was found in previous studies 10,[20][21][22][23]47 . An attempt at the kinetic contribution (diffusion coefficient) failed due to force field limitations. Therefore, this study cannot explain the low temperature anomaly in the work of formation, except that the thermodynamic part has been confirmed. This work does illustrate the possible competing pathways that are available during a nucleation event. The work of formation is dependent on the nucleating pathways which makes Eq. (1) inherently more complex. As a consequence, Volmer and Weber's kinetic model 30 would have different activation energies for adjoining clusters based on their thermodynamic reactant state shown here. Even for a simple binary glass, nucleation is very complex and additional studies are needed. Once scientific investigation answers a question, it unlocks many more questions. For instance, does silica clustering drive glass-ceramic nucleation? Do all polymorphs of barium silicate exist in other compositions? What does the core tell us about nucleation? How does the silica rich core develop as a function of time and temperature? Can the core's creation be perturbed by adding a fast diffusing alkali metal to the system? Can a phase sensitive glass force field be generated? Further work into advancing Diffuse Interface Theory and Density Functional Theory methods, such as competing pathways, regions of slow mobility, and interface development with advancements in experimental methodologies, will increase our understanding of nucleation.

Computational methods
The results presented here are based on studies from a Grand Canonical Monte Carlo (GCMC) model, molecular dynamics (MD) models, and experimental observations. The GCMC model is used to explore the nucleation free energy (NFE) landscape of cluster sizes and compositions, a technique that has shown to give good agreement with experimental results 26 . MD is used to probe the crystal-melt interface, which is an atomistic analogue to DIT. One item the GCMC model does not account for is diffusion, so the MD model is used to investigate this contribution. Experimentally, the nucleated barium silicate glasses were analyzed using Scanning Electron Microscopy (SEM), Transmission Electron Microscopy (TEM), and Energy Dispersive Spectroscopy (EDS). This section outlines each of these models, analysis techniques, and the experimental methods used. 26,48 was developed to study nucleation in a lithium disilicate glass. The work of cluster formation is not described by the usual surface tension and thermodynamic driving terms but instead tracks (1) the amount of work that is required to form a cluster, (2) the energy required to crystallize that cluster and (3) the associated solvation energy. The model only simulates the cluster; it does not simulate the atoms in the glass. This allows the computation to focus solely on the cluster development instead of using computational resources to move atoms in the glass region that do not directly partake in the nucleation event. Each of the clusters experiences the continuum field of the glass in the form of an implicit solvation model. For completeness some details of the simulation technique are discussed below. Full details of the computer code and methods can be found in Refs. [49][50][51][52][53] .

Grand canonical Monte Carlo (GCMC) model. This model
Cluster formation. This model uses the Grand Canonical (μVT) ensemble, where the clusters can have a fluctuating number of atoms. The Grand Canonical ensemble has the cluster physically separated from (but thermodynamically coupled to) the liquid phase. This is similar to the Gibbs approach except it does not explicitly model the liquid phase. We only track the development and growth of one cluster in the Grand Canonical ensemble. It is difficult to directly obtain the chemical potential, μ, for these glass-ceramic systems, whether one uses a model approach (e.g., Widom particle insertion technique 54 ) or experiment. Instead, from statistical mechanics, one can derive the chemical potential to be equivalent to the product of the Boltzmann constant, absolute temperature, and the natural logarithm of the number density 55 . This constant volume term in the ensemble is enforced by an energy-based Stillinger-type cluster criterion. The cluster is defined as a group of atoms in which every atom has at least one neighbor with interaction energy equal to or greater than 50% of the Si-O bond energy (226 kJ/mol). The atomistic models use the force field developed by Pedone et al. 56 . All clusters from 1-40 SiO 2 units and 0-25 BaO units have been sampled with at least 800 million Monte Carlo cycles for temperatures of 980 K, 1090 K, and 1160 K (representing below, at, and above T max ).
We utilize the Aggregation-Volume-Bias Monte Carlo (AVBMC) move type 49 to bypass the low acceptance rate of a traditional Metropolis Monte Carlo move 35 in these simulations. An AVBMC move has two options: (1) cluster addition (out → in), or (2) cluster subtraction (in → out). The target particle j is then randomly selected in the cluster. If this is a cluster addition move, the code will randomly select a particle from the ideal glass phase and place it inside the V in volume near particle j. In the case of a subtraction move, a random particle within the V in region of particle j is chosen as a candidate for deletion. The energy change between the current configuration A and the trial configuration B (based on addition or subtraction) is then calculated, i.e., ΔE = E B − E A . Within the V in region of the cluster N in is the total number of particles (sums to N − 1, where N is the total particles of the cluster, as j is excluded). The trial move is accepted with the following probabilities: 1. Cluster addition 2. Cluster destruction www.nature.com/scientificreports/ Incorporating the energetic and entropic factors experienced in the associating and non-associating volumes will lead to a higher probability for accepting a trial move, hence increasing the efficiency of the Monte Carlo moves and leading to an effective evolution of the cluster distributions from pre-to post-critical sizes. The Monte Carlo move types are chosen at a frequency of 30% AVBMC and 70% translational movements. The resulting cluster structure and calculated energies are invariant on the Monte Carlo move frequencies, however the total simulation length does depend on the move ratio. It was found that an increased amount of translational movement tends to help the cluster re-organize whenever a particle is added/removed. Crystal transition. Steinhardt order parameters 9,57 were used to measure the amount of energy required to move the atoms in the amorphous structure into the crystalline atomic positions. This technique has been used for liquid to solid and solid to solid phase transitions [58][59][60] . At the heart of this computational technique, the bond vectors and angles are calculated and compared to reference crystal values. The differences in the order parameter values can be sampled to calculate the amount of energy required to form a crystal from the initial amorphous positions 26 .
Solvation energy. The solvation energy is found using the Implicit Glass Model (IGM) 61 . IGM is a member of the Generalized Born implicit solvation model 62,63 family that are normally used for proteins. This is a per atom solvation term that is incorporated into the cluster's total energy. The parameters used in this work came from Ref. 61 . While this term typically only adds ~ 1 k B T unit to the work of cluster formation, it is very important since the sampling of the cluster formation is conducted only on one cluster in vacuum (and no residual glass).
Nucleation free energy (NFE) analysis. AVBMC's purpose is to increase the frequency of cluster growth and destruction event. This is still insufficient when attempting to sample clusters near the critical size. Due to the low occurrence probability of clusters at the critical size, a typical barrier is of the order 30 k B T; this means that one must sample e 30 times to observe the event! This is solved by using an umbrella sampling scheme 64 in which a bias potential is developed during the simulation to enhance the sampling frequencies. The bias potential is removed in the nucleation free energy analysis.
In the Grand Canonical ensemble, the probability, P, of observing a particular state described by n particles and their configuration r n in a biased distribution is where β is 1/k B T (where k B is the Boltzmann constant and T is the temperature in absolute units and B(n) is the biasing potential as a function of cluster size n. The bias potential is multiplied by the AVBMC acceptance criteria (Eq. 4, 5) as W(n + 1)/W(n) (for addition) or W(n − 1)/W(n) (for cluster destruction). The probability of having clusters of size n in the bias distribution, P(n) biased is computed; the probability of the unbiased distribution is then given by The objective is to choose a proper weighting function, to have a statistically precise estimate of P(n) biased with the least amount of computational time. A reasonable criterion is to have an equal probability of observing each cluster size n. This can also provide the maximum possible sampling of each size for a given simulation length.

The choice of this weighing function is
This gives an equal probability to sample all cluster sizes of n in the biased distribution. The free energy of cluster formation with this biasing potential becomes Unfortunately, this choice generates a new problem, as the unbiased probabilities are not known (and their determination is the objective of the simulation). There are different ways to obtain the appropriate biasing potential. From knowledge of this algorithm, there are two reasonable starting choices: (1) conduct a few million samplings of unbiased small cluster sizes to initiate the bias, or (2) if another bias potential is known (either generated at another temperature or extrapolated from the previous point) those values can be used. Since this bias approach is an iterative and adaptive process, many millions of Monte Carlo cycles are needed to ensure convergence. Convergence is reached when the bias potential does not change by more than 1% of its value. After every 10,000 Monte Carlo cycles the bias potential is re-calculated. In total, each cluster size requires only   (6) P n, r n = exp (µnβ) exp −βE N, r n B(n) , (7) P(n) unbiased = P(n) biased W(n) . www.nature.com/scientificreports/ one CPU core for approximately 2-3 weeks of clock time. It should be mentioned that for the lowest temperature studied (980 K), this convergence required 30% more Monte Carlo cycles than for the other temperatures.
Molecular dynamics. The Monte Carlo model allows one to probe the complete nucleation pathway. A molecular dynamics model is used to study the explicit interface between the crystal and the glass, the evolving melt structure near the interface, and the bulk diffusion of the barium silicate glass. A stoichiometric barium disilicate glass was first generated by randomly inserting 250,000 atoms into a cubic simulation box and equilibrating them at 3000 K for 2 ns in the Canonical (NVT) ensemble. This liquid was then continuously cooled down to 300 K over a period of 6 ns in the Isothermal-isobaric (NPT) ensemble under atmosphere pressure. The final sample dimension was 23.2 nm × 46.3 nm × 2.89 nm. Afterwards, a spherical void with a radius of 4.1 nm was created in the center of the glass by deleting atoms inside the region and refilling it with a piece of pre-made spherical BaO 2SiO 2 crystal of the same radius. The premade crystal was cut out of a large piece of crystal generated using the crystal structure information reported in Refs. 65,66 . Because of the atom deletion and insertion operations, the system usually lost charge balance. To remedy this, additional atoms in the glass were randomly chosen and deleted to maintain charge neutrality for the final sample. The bulk diffusion study used the same sized simulation box and heat treatment, but without the inserted crystal.
To monitor the evolution of the glass-crystal interface, the sample was gradually heated up from 300 K to the target temperature (960 K, 1040 K, 1190 K) in the period of 1 ns using the NPT ensemble, and further held at the target temperature for 1 ns. The bond-orientational order parameter Q4, introduced by Steinhardt et al. 57 , was calculated to characterize the local orientational order in the atomic structure. This parameter was expected to capture the structural differences between the glass, the interface, and the crystal. To obtain a better statistical Q4 value, the atoms were spatially binned into spherical shells from the center of the simulation box to the edge, using a bin size of 0.3 nm.
These simulations were conducted using LAMMPS 67 with a timestep of 2 fs. Temperature and pressure were well controlled via the Nose-Hoover 68,69 thermostat and barostat, respectively. Periodic boundary conditions (PBC) were applied in all directions over the course of these operations. The force field created by Pedone et al. 56 was used to capture the short-range interactions for both phases. For the long-range columbic interactions, the damped shifted force (dsf) method, with a cutoff of 8 Å and damping parameter of 0.25 Å −1 , was used to speed up calculations.

Experimental methods
Glass ceramic samples were prepared by Corning Incorporated using the melting and quenching procedures discussed by Xia et al. 21 . The source materials were barium carbonate and silica.
TEM images were acquired using a bright-field detector in an aberration corrected FEI Titan high resolution transmission electron microscope (HRTEM) operated at 200 keV electron beam energy in the scanning mode (STEM). The samples were prepared by the conventional TEM preparation method (mechanical polishing and dimpling followed by Ar-ion polishing).
Samples were prepared for SEM analyses by fracturing and acid etching in a 0.1% HF solution for 10 s to help reveal surface topography. An amorphous conductive carbon coating was evaporated onto the samples to reduce charging. Secondary electron images were acquired in a Zeiss 1550 field emission scanning electron microscope at an accelerating potential of 5 keV and 200 picoamps beam current. Energy-dispersive X-ray spectroscopy analyses were acquired using an Oxford Instruments 80 mm 2 silicon drift detector and AZtec X-ray microanalysis software.