{'fraction_non_alphanumeric': 0.03584964706700511, 'fraction_numerical': 0.01436072186098265, 'mean_word_length': 4.692794932699921, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polymorph selection during crystallization of a model colloidal fluid with a free energy landscape containing a metastable solid

The free energy landscape responsible for crystallization can be complex even for relatively simple systems like hard sphere and charged stabilized colloids. In this work, using hard-core repulsive Yukawa model, which is known to show complex phase behavior consisting of fluid, FCC and BCC phases, we studied the interplay between the free energy landscape and polymorph selection during crystallization. When the stability of the BCC phase with respect to the fluid phase is gradually increased by changing the temperature and pressure at a fixed fluid-FCC stability, the final phase formed by crystallization is found to undergo a switch from the FCC to the BCC phase, even though FCC remains thermodynamically the most stable phase. We further show that the nature of local bond-orientational order parameter fluctuations in the metastable fluid phase as well as the composition of the critical cluster depend delicately on the free energy landscape, and play a decisive role in the polymorph selection during crystallization.


I. INTRODUCTION
The remarkable argument of Alexander and McTague, based on the Landau theory, that in the case of a simple fluid undergoing weak first order phase transition, the body-centered-cubic (BCC) phase should be favored over the thermodynamically stable face-centered-cubic (FCC) phase [1] has till now defied a quantitative understanding. It has been hard to reconcile with the classical nucleation theory (CNT) [2][3][4][5], which in its simplest form, considers the competition between the free energy gain in the fluid to solid transformation and the cost of creating the fluid-solid interface; the net free energy change controls the crystallization process. CNT may capture the kinetics for the cases where only two free energy basins (the metastable parent and the stable daughter phases) are involved, but it could fail to address phase transformation in the presence of multiple solid phases.
The participation of the intermediate metastable phase(s) through wetting of the stable phase nuclei is an important factor that needs to be taken into account. In such cases, we have a competition between thermodynamic (the stability of solids) versus kinetic (the free energy barrier of nucleation) control. Empirically, one can invoke the Ostwald's step rule [6] that states that the crystal phase that forms out of metastable melt need not be the thermodynamically most stable phase, rather it is the one closest in stability to the parent phase. However, the Ostwald's step rule still lacks a solid theoretical foundation. Stranski and Totomanow [7] argued that the solid with lowest free energy barrier will emerge from the metastable fluid irrespective of its stability with respect to other solid phases, which is yet to be demonstrated in a fully microscopic computational study.
Recent advancement in experimental techniques has led tremendous interest in understanding and controlling the phase transformation and polymorph selection in complex materials and biological systems [8][9][10][11][12][13][14][15][16][17]. Complex materials are often characterized by the presence of multiple length and energy scales in the interaction potential between the constituents, and usually display complex energy landscape [18] and rich phase behavior consisting of multiple phases [19][20][21][22][23][24][25][26]. Understanding the nature of the free energy landscape and its connection with the pathways of phase transition has great practical relevance in seemingly distinct branches of science ranging from materials (e.g. polymorph selection [8,27]) to biological (e.g. protein crystallization, aggregation [12,28,29]) sciences. Crystal engineering relies heavily on gaining a molecular level understanding and control of the free energy landscape, and in turn, the pathways of phase transformation [12,17]. Multiple intermediate metastable phases are already known to play an important role in the formation of random spin and structural glasses in both experimental and theoretical descriptions [30][31][32][33].
In experiments, proper characterization and controlled change of the complex free energy landscape pose a major challenge in gaining fundamental understanding of the crystallization and polymorph selection processes. For example, any change in thermodynamic conditions such as temperature (T ) and pressure (P ) alters the whole free energy landscape (the stability of all the phases), and thus prevents us from understanding how the relative stability of a particular intermediate metastable polymorph would affect the pathways and rate of crystallization. Using phenomenological classical density functional theory (DFT) [5], we recently showed that the controlled changes of the (meta)stability of intermediate phases can give rise to diverse non-classical pathways of phase transformation, ranging from wetting-mediated to Ostwald's step rule-like scenario [27,34]. However, precise microscopic pathways and criteria for polymorph selection in such complex systems largely remain elusive and demand further controlled studies with atomic resolutions.
In this work, we used hard-core repulsive Yukawa as a model system to computationally explore the polymorph selection during crystallization. This system exhibits a rich phase behavior consisting of fluid-FCC/BCC and FCC-BCC phase coexistence lines along with two triple points [35], and thus, ideally suited to study the interplay between the free energy landscape and the selection of FCC and BCC polymorphs. We explored the diverse pathways of phase transformation through the controlled change of the free energy landscape. On gradually increasing the stability of the BCC phase with respect to the fluid phase at a fixed fluid-FCC stability, we observed a cross-over from the formation of the thermodynamically most stable FCC phase through a wettingmediated pathway to an Ostwald's step rule-like scenario where the BCC phase of intermediate stability (stable with respect to the fluid and metastable with respect to the FCC) grows despite FCC being the thermodynamically the most stable phase. We further observed that the composition of the critical cluster depends delicately on the free energy landscape, and plays a decisive role in the polymorph selection during crystallization. Additionally, we also explored the microscopic pathways of the emergence of the composition of different (FCC and BCC) polymorphs in the critical cluster from the metastable fluid phase.

A. Model details
We performed Monte-Carlo (MC) simulations [36] on a system interacting via hard-core repulsive Yukawa potential, where σ is the particle diameter and is the energy at contact distance σ. β = 1/k B T where k B is the Boltzmann constant and T is the temperature. σ and are used as units of length and energy, respectively. In this study, we have truncated the interaction potential at a distance r c = 3.0σ and shifted to zero, and chose kσ = 5 as this value is commonly used in majority of the studies on Yukawa system [14,35], and the phase diagram for this parameter choice suits best to this work (note that, the phase behavior of Yukawa system shows strong dependence on the choice of the potential parameters [37]).

B. Computation of phase diagram
To obtain the phase diagram, we first computed the Helmholtz free energy (F ) of FCC and BCC phases at β = 8 and ρσ 3 = 0.75 using thermodynamic integration [36] in NVT ensemble for a system consisting of N = 432 particles for BCC and N = 500 particles for FCC. Using these Helmholtz free energies, we obtained the chemical potential at different pressures by integrating the equation of state [36]. The chemical potential of the fluid phase was computed using Widom insertion method [36] at β = 8 and ρσ 3 = 0.15. Again, integration over the equation of state was performed to determine the chemical potential of the fluid phase as a function of pressure. Equating these chemical potentials of fluid and solid (BCC and FCC) phases we obtained the coexistence pressures βP σ 3 = 26.5 and 27.4 at β = 8 for fluid-BCC and fluid-FCC, respectively. We verified these coexistence pressures using direct free energy calculation employing umbrella sampling [38] and obtained βP σ 3 = 26.8 and 27.7 for fluid-BCC and fluid-FCC, respectively. After computing the coexistence fluid-solid pressures at β = 8, we obtained the fluid-BCC and fluid-FCC coexistence pressures as a function of β using the Clausius-Clapeyron relation [39,40], where ∆h = h j (β 1 , P 1 ) − h i (β 1 , P 1 ) is the enthalpy difference per particle between phases j and i, and ∆v = v j (β 1 , P 1 ) − v i (β 1 , P 1 ) is their volume difference per particle at inverse temperature β 1 and pressure P 1 . P 1 is the coexistence pressure at inverse temperature β 1 , and P 2 is the coexistence pressure at inverse temperature β 2 . We carried out NPT MC simulations with 20, 000 equilibrium MC steps (1 MC step equals N numbers of single particle displacement and one volume move attempts) followed by 50, 000 production steps at β 1 and P 1 in both the phases, i and j, and computed ∆h and ∆v. The systems contained 432 particles for fluid-BCC and 500 for fluid-FCC coexistence lines. Using Eq. 2, we obtained fluid-BCC and fluid-FCC coexistence pressures as a function of β and they cross each other at a triple point, β = 4.5 and βP σ 3 = 26.2. Starting from this triple point we computed the BCC-FCC coexistence line as a function of β using Clausius-Clapeyron equation as described above. Finally, combining these three coexistence lines -fluid-BCC, fluid-FCC and BCC-FCC -we obtained the phase diagram shown in Fig. 1. The chemical potential differences of FCC and BCC phases with respect to the fluid phase at thermodynamic conditions studied in this work (asterisks in Fig. 1, and Table I) were computed using thermodynamic integration for larger systems consisting of N = 2662 and 2916 particles for BCC and FCC phases, respectively.

C. Identification of solid-like particles and polymorphs
Solid-like crystallites in the metastable fluid phase were identified using the method introduced by Frenkel and co-workers [41]. This method first identifies the local bond-orientational symmetry of particle i using a complex vector q lm (i) [42] as , where N b (i) is the number of nearest neighbors of the i th particle. Two particles were considered to be neighbors if the distance between them (|r ij |) was less than the cut-off distance of q c = 1.38/(ρσ 3 ) 1/3 , where q c is the radius of the first shell of FCC lattice (measured from the position of the minimum separating the first and second peaks in the radial distribution function), ρσ 3 is the reduced density. Y lm (r ij ) is the spherical harmonics and r ij is the distance vector between the particle i and its neighbor j. l and m are integers with −l ≤ m ≤ l. The unit vector of q lm (i) is given by, Using the unit vector d lm (i), a scalar product S l (i, j) which measures the correlation in bond orientational order between neighboring particles can be defined as, where the superscript * indicates complex conjugate. Two neighboring particles i and j are considered to be connected if S 6 (i, j) > 0.7. The particle i is identified as solid-like if the number of such connections is more than 7.
In order to assign the polymorphic identity of a solidlike particle, we employed the coarse-grained (or locally averaged) bond orientational order parameter introduced by Lechner and Dellago [43]. Using the order parameter given in Eq. 3, one can define a locally averaged complex vectorq lm (i) as, where j = 0 indicates the particle i itself. Given the coarse-grained complex vectorq lm (i), one can further define coarse-grained order parametersq l (i) andw l (i) as and where the term in the parentheses (...) indicates the Wigner 3j symbol. The integers m 1 , m 2 and m 3 range from −l to l and only the terms with m 1 + m 2 + m 3 = 0 are allowed to contribute to the summation. Oncew l (i) is defined, we identify a previously assigned solid-like particle as BCC-like ifw 6 > 0, whereas it is considered to be HCP-like ifw 6 ≤ 0 andw 4 > 0. A particle is considered as FCC-like ifw 6 ≤ 0 andw 4 ≤ 0 [44].

D. Computation of nucleation free energy
The fluid to solid nucleation free energy profiles shown in Fig. 2 were computed by employing umbrella sampling method [45] with size of the largest cluster as the order parameter in NPT ensemble consisting of N = 2916 particles. The force constant of the umbrella potential was taken to be λ = 0.1k B T . The fluid to FCC nucleation free energy profile was computed by biasing the system along the size of the largest FCC cluster with force constant λ = 0.1k B T , while simultaneously preventing the formation of BCC clusters using another umbrella potential with λ = 100k B T along the size of the largest BCC-like cluster having minimum at cluster size n b = 0. Similarly, for fluid to BCC nucleation free energy profile, two harmonic biasing potentials -one along the size of the largest BCC cluster with λ = 0.1k B T and the other with λ = 100k B T along the size of the largest FCC cluster having minimum at cluster size n f = 0 -were employed simultaneously.

A. Phase diagram
In Fig. 1, we show the computed phase diagram of the kσ = 5 hard-core repulsive Yukawa system (for details of the phase diagram computation see Section II B). The red asterisks in Fig. 1 denote the thermodynamic conditions at which crystallization has been studied in this work. At these conditions, the fluid is metastable with respect to both BCC and FCC phases. FCC phase is the thermodynamically most stable and the BCC phase is stable with respect to the fluid phase, but metastable with respect to the FCC phase (Table I). As one moves from I to III via II in the phase diagram (Fig. 1), the free energy difference between the fluid and FCC remains same, however, BCC gradually becomes more and more stable with respect to the fluid phase (Table I) Table I, at these state points, the chemical potential difference between the fluid and the FCC solid (β∆µFCC) is fixed at −0.28 and the stability of BCC phase with respect to the fluid phase (|β∆µBCC|) gradually increases on moving from I to III. To uncover the role of the presence of metastable BCC phase of varying stability on the microscopic mechanism of crystallization, in Fig. 2 (top panel), we show the dependence of the composition of the largest solid-like cluster (fraction of FCC, BCC and HCP-like particles in the largest solid-like cluster) on its size. We follow the method introduced by Frenkel and coworkers [46] to define solid-like particles based on local bond-orientational order parameter q 6 [42]. The solid-like particles in the largest cluster are further identified as BCC, FCC and HCP-like based on their averaged localw 6 andw 4 order parameters [43,44] (see Section II C for the details). A solid-like particle withw 6 > 0 is considered as BCClike, whereas a particle withw 6 ≤ 0 andw 4 > 0 as HCP-like. FCC-like particles are those withw 6 ≤ 0 and w 4 ≤ 0. This useful assignment criteria to distinguish different solid-like local environments in the metastable fluid phase has been used extensively in computer simulations of model atomic and molecular systems [47][48][49][50].
At the thermodynamic conditions studied here, the metastable fluid does not undergo spontaneous phase transition on the simulation time scale. Therefore, we employed metadynamics simulations [51] considering the size of the largest solid-like cluster as order parameter to assist the system to overcome the nucleation free energy barrier and grow spontaneously [52,53] (note that, as the bias is on the size of the largest solid-like cluster, it does not affect the natural selection of polymorphs). When BCC is only marginally stable with respect to the fluid phase (at I), we observe nucleation of the FCC-like clusters. On increasing the stability of the BCC phase (at II), we observe a competitive growth of the both, FCC and BCC-like clusters. On further increasing the stability of the BCC phase (at III), we observe an Ostwald's step rule like scenario where BCC phase of intermediate stability nucleates from the fluid, despite FCC being the thermodynamically most stable phase (Fig. 2, top panel). This observed crossover from the FCC-dominated to the BCC-dominated cluster on gradual increase of the stability of the BCC phase is consistent with the predictions of our recent classical DFT [27]. As the extent of HCP-like particles at all the three conditions -I, II and III -is low, from now onwards FCC denotes FCC+HCP unless HCP is explicitly specified.
Furthermore, by using umbrella sampling [36,45], we also computed the nucleation free energy barrier for the solid-like clusters along with the FCC and BCC-like clusters (Fig. 2, bottom panel). For the FCC nucleation free energy profile, we first identified FCC and BCC-like particles using the above mentioned criteria and then suppressed BCC-like fluctuations by imposing an umbrella bias along the number of BCC-like particles (n b ) with its minimum at n b = 0. For the BCC nucleation free energy profile, the same approach was followed except that the FCC-like fluctuations were suppressed, in place of BCC (see Section II D for the details). This approach provides conditional free energy barrier of nucleation for the pure FCC and BCC phases where the (indirect) participation of the other phases -BCC and FCC, respectively -is either absent or negligibly small. This conditional free energy barrier can provide an estimate of the change in the free energy barrier due to compositional heterogeneity of the critical cluster.
As evident from Fig. 2, the free energy barrier of crys- tallization decreases on moving from I to III, even though the stability of FCC (β∆µ FCC ) remains the same. The fluid-to-FCC nucleation barrier (in absence of wetting by BCC-like particles) remains approximately the sameconsistent with the fixed β∆µ FCC at all the three conditions (Table I). Due to the gradual increase of the stability of the BCC phase on moving from I to III, the fluid-to-BCC nucleation free energy barrier deceases and crosses the fluid-to-FCC nucleation barrier near the thermodynamic condition of II. This crossover in the nucleation barrier leads to a transition in the nucleation mechanism from the formation of thermodynamically favored FCC phase to the formation of the metastable BCC phase (Fig. 2, top panel), unambiguously suggesting an Ostwald's step rule like scenario.
The polymorphic identity of the final solid formed after crystallization from the metastable fluid phase at state points I, II and III (Fig. 2) was confirmed by comparing the radial distribution functions (RDFs) and the local bond-orientational order parameters (w 4 andw 6 ) of the solid phases with respective pure FCC and BCC phases equilibrated at the same thermodynamic conditions (shown in the top and bottom panels of Fig. 3, respectively). At I and III, the RDFs and the localw 4 the top panel and red circles in the bottom panel) phases, respectively. At II, however, the final solid phase can be considered as a mixture of FCC and BCC phases, which is a consequence of competitive nucleation and growth at this thermodynamic condition (Fig. 2).

C. Composition of critical clusters
To further uncover how the composition of the critical cluster depends on the free energy landscape or the relative stability of different solid phases, in Fig. 4, we show representative snapshots along with composition profiles of the critical cluster averaged over hundred independent simulations at I, II and III. As the snapshots show, at I, the critical cluster is composed of mostly FCC-like particles with random patches of BCC-like particles at the surface only partly covering the nucleating and growing cluster. At II, the critical cluster is a mixture of both FCC and BCC-like particles of comparable fraction, and at III, we observe BCC-dominated cluster along with dispersed FCC and HCP-like particles.
At the bottom panel of Fig. 4, we show the composition profiles, defined as the number of particles of ith solid phase in a shell of radii r and r + ∆r divided by the total number of solid-like particles in that shell, for FCC, BCC and HCP-like particles in the critical clusters along with the normalized density of the solid-like  I I I I I   0  2  4  6   II II II II II   0  2  4  6   III III III III  particles (black line) as a function of the distance from the center of the cluster. At I, the fraction of BCC-like particles increases on moving outward from the center to the surface (denoted by the vertical dashed line) of the cluster, indicating preferential wetting of the surface by the BCC phase. At II, the FCC and BCC-like particles are randomly distributed throughout the cluster in a similar proportion. At III, on the other hand, BCClike particles dominate throughout the cluster and we do not observe any signature of wetting of the surface by the FCC phase. The lower free energy barrier of crystallization compared to the free energy barriers of pure phases (Fig. 2, bottom panel) can be attributed to the compositional heterogeneity and wetting of the solid-like clusters.
Although recent computer simulation studies on atomic systems (such as Lennard-Jones, Gaussian core model, hard-core Yukawa and hard sphere) [44,47,48,50,[54][55][56][57] demonstrate various pathways of crystallization at different thermodynamic conditions, an exact criterion for the change in mechanism of crystallization from wetting-mediated to Ostwald's step rule was still missing. In the present study, controlled change of the free energy surface (or more precisely, on changing the stability of the BCC phase at a fixed stability of the FCC phase with respect to the fluid) and the explicit computation of the conditional free energy of the pure polymorphs (Fig.   2) enabled us to gain quantitative understanding of this criterion.

D. Fluctuations of local bond-orientational order parameter in the metastable fluid
In a metastable fluid, through thermal fluctuations, crystallites of relatively stable phase(s) appear and disappear and sometimes grow leading to the phase transition. Recent studies show that the key to polymorph selection is hidden in the bond-orientational order parameter fluctuations in the metastable fluid [44,48,58]. To explore how the stability of the intermediate BCC phase alters the local structural fluctuations in the metastable fluid phase, in Fig. 5, we show the computed localw 6 distribution (P (w 6 )) along withw 4 distribution of the particles withw 6 ≤ 0 at metastable state conditions I, II and III. In computation of P (w 6 ) and P (w 4 ), we consider particles withq 6 > 0.27 [44]. The former enables us to distinguish BCC and HCP/FCC-like, and the latter distinguishes HCP and FCC-like local structural fluctuations. On gradually increasing the stability of the BCC phase (I→III), we observe suppression of FCC/HCP-like fluctuations and concurrent enhancement of BCC-like fluctuations (Fig. 5). In Fig. 6 (top panel), we show the evolution of the localw 6 distribution with the extent of crystallinity which is quantified using the localq 6 order parameter in the metastable fluid phase. At I, the system becomes increasingly enriched with FCC/HCP-like environments (w 6 < 0) on including the particles with the higherq 6 values. A similar behavior is observed at II, which is consistent with the recent observation for the Gaussian Core Model system [48]. In contrary, at III, the population is slightly biased towards BCC-like environments (w 6 > 0). To further analyze the polymorphic identity of the particles withw 6 < 0, in the bottom panel of Fig. 6, we present the distribution of localw 4 for the particles withw 6 < 0 on varying the range ofq 6 (same as in the top panel). This order parameter distinguishes the FCClike (w 4 < 0) and HCP-like (w 4 > 0) local environments. As the figure suggests, at all the three conditions, the distributions are biased towards FCC-like environments (w 4 < 0). That is, the particles withw 6 < 0 are dominated by the FCC-like local environments (w 4 < 0), as is also evident from Figs. 2, 4 and 5(b). The results shown in Fig. 5 and Fig. 6 highlight the underlying connection between the free energy cascade and the nature of the fluctuations in the metastable fluid and confirm recent observations that the information about the polymorph selection is indeed encoded in thermal fluctuations of local bond-orientational order parameters in the metastable fluid phase. Finally, as the composition of the clusters controls the free energy cost of their formation ( Fig. 2 and Fig. 4), we address the question of the selection and emergence of the critical cluster composition (fraction of different solid-like particles in the critical cluster) by fluctuations off the metastable fluid phase. Following the history of each solid-like particle in the critical cluster one can gain a mechanistic understanding of how each of the solidlike particles in the critical cluster eventually form from fluid-like environment. The transformation of a fluid-like . Note that at I, FCC-direct path is the dominant contributor, and at III, BCC-direct path is the dominant contributor to the formation of solid-like particles in the critical cluster. At II, we observe competitive direct appearance of FCC and BCC-like particles from respective fluid-like environments. The transition fractions were obtained by averaging over 100 independent trajectories. The lower panel depicts a schematic representation where grey, blue and red spheres indicate fluid, FCC and BCC-like particles, respectively, and the width of the arrow for each step is proportional to its weight.
particle to the solid-like can occur either through (a) direct transformation of the fluid to the final solid-like environment or (b) indirect pathway where fluid-like particle first transforms to BCC(FCC) and then to FCC(BCC).
We computed the weights of these two pathways for each of the solid-like particle belonging to the critical cluster (Fig. 7). At I, the majority of the solid-like particles in the critical cluster (∼ 60%) forms via the direct fluid→FCC pathway and only 20% via the indirect (fluid→BCC→FCC) pathway (Fig. 7, top panel). On the other hand, at III, the majority of transitions (∼ 50%) is via the direct fluid→BCC pathway. At II, however, we observe a competitive direct appearance of FCC and BCC-like particles (ca. 30% and 40%, respectively). In all the three cases, the solid-like particles formed via indirect pathways are only 15-20%. This observation suggests that the composition of the critical cluster is predominantly guided by the direct transformation of the fluid-like particles to the solid-like rather than indirect pathways mediated by solid-solid transformations. The lower panel of Fig. 7 depicts a schematic representation of the whole process.

IV. CONCLUSIONS
Using hard-core repulsive Yukawa as a model chargestabilized colloidal system, we demonstrate here the kinetic origin of preferential formation of the BCC phase even though FCC is thermodynamically the most stable phase, thus unambiguously justifying the age-old observation of Alexander and McTague [1]. In this process, this study brings out the true essence of the Ostwald's step rule [6,7]. We further show that the nature of the local bond-orientational oder parameter fluctuations in the metastable fluid phase as well as the composition and size of the critical nucleus depend delicately on the relative stability of the intermediate BCC phase. The composition of the critical cluster is guided by the direct transformation of fluid-like particles to the solid-like rather than indirect pathways involving solid-solid transformations. In addition, the results obtained in this work qualitatively validate the predictions of our recent phe-nomenological classical DFT [27]. As this theory is not specific to any particular system or inter-particle interaction potential, we anticipate that, irrespective of the system, if the nature of the free energy landscape is like the one considered here, the results obtained in this work should hold true. We also anticipate that this study could provide important insights into the synthesis of polymorphs of desired structures and properties by controlled change of the free energy landscape either through changing the thermodynamic conditions or altering the interparticle interactions (especially in macro-and mesoscale systems) [26,[59][60][61].