An investigation of the kinetics and thermodynamics of NaCl nucleation through composite clusters

Abstract Having a good understanding of nucleation is critical for the control of many important processes, such as polymorph selection during crystallization. However, a complete picture of the molecular-level mechanisms of nucleation remains elusive. In this work, we take an in-depth look at the NaCl homogeneous nucleation mechanism through thermodynamics. Distinguished from the classical nucleation theory, we calculate the free energy of nucleation as a function of two nucleus size coordinates: crystalline and amorphous cluster sizes. The free energy surface reveals a thermodynamic preference for a nonclassical mechanism of nucleation through a composite cluster, where the crystalline nucleus is surrounded by an amorphous layer. The thickness of the amorphous layer increases with an increase in supersaturation. The computed free energy landscape agrees well with the composite cluster-free energy model, through which phase specific thermodynamic properties are evaluated. As the supersaturation increases, there is a change in stability of the amorphous phase relative to the solution phase, resulting in a change from one-step to two-step mechanism, seen clearly from the free energy profile along the minimum free energy path crossing the transition curve. By obtaining phase-specific diffusion coefficients, we construct the full mesoscopic model and present a clear roadmap for NaCl nucleation.


Introduction
Crystallization from solution is an important process widely used in the pharmaceutical, chemical, and food industries for separation and purification. Controlling these processes to achieve the desired form of the end product requires a good understanding of the underlying mechanism of nucleation. However, our understanding of nucleation is still far from complete, as a result of the experimental challenges for observing nucleation due to short time and length scales involved in the formation of a critically sized nucleus except for a few systems that are slow enough to be captured by simple microscopy (1,2).
Classical nucleation theory (CNT) offers a simple phenomenological description of nucleation, and has been used widely in theoretical (3)(4)(5) and experimental (6)(7)(8) studies of nucleation. CNT is based on the assumption that nucleation occurs via a single-step process, where the system overcomes a single free energy barrier and the clusters forming in the solution have the same properties as that of a bulk crystal (9,10). Thus, the size of the emerging nucleus, n, is the only slow variable and enough as a reaction coordinate to describe the whole nucleation process. However, nucleation rate predictions of CNT can deviate from experiments by several orders of magnitude (11). Nonclassical nucleation mechanisms have been pointed out as one possible explanation for the discrepancies. In fact, there is evidence for various nonclassical mechanisms like multistep mechanism (12), and nucleation through prenucleation clusters (13,14), composite clusters (15,16), or polycrystalline structures (17). The common feature of these nonclassical mechanisms is that the cluster goes through some structural change as it is emerging from the mother phase and growing to reach the critical size. Thus, some structural order parameter is required in addition to the nucleus size to describe such transformations.
While nucleation is difficult to access experimentally, simulations can provide valuable atomic level insight into the mechanism of nucleation from solution. NaCl nucleation from aqueous solution has been the focus of many computational studies (5,(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29) since it is a simple, well-parameterized system with force fields for which the predicted solubility values are known (30,31). Conclusions on the nucleation mechanism for this system have been varied: large-scale MD simulations done by Chakraborty and Patey (21,22) showed evidence for a two-step mechanism, where nucleation originated in a region of higher salt concentration than that of the bulk solution. They reported early stage nuclei consisting of loosely ordered ions with significant amount of waterinconsistent with the capillarity assumption of CNT. A later study that also employed direct simulations of NaCl nucleation (23) revealed that the expected lifetime and nucleation probability of clusters depended not only on their size but also the geometric arrangement of their ions. Zimmerman et al. used seeded simulations to estimate the nucleation rates at 3 different supersaturations and were able to match the trend of experimental measurements. It is worth noting that the seeding approach uses elements from CNT for calculating the nucleation rate, inherently assuming a single-step nucleation mechanism (32). A later study (25) showed that the chemical potential of NaCl ions in solution increased as a function of concentration, until it reached a plateau after a concentration of 15.0 mol/kg.This plateau is explained with a spinodal decomposition that leads to a barrierless, spontaneous formation of amorphous clusters, which is the first step of the nucleation process. Then, crystallization within the amorphous clusters proceeds via a finite barrier, which constitutes the second step.
In this study, we examine the mechanism of NaCl nucleation from solution and provide a complete thermodynamic picture that reveals the preferred pathway of nucleation. As pointed out by earlier studies (23), the degree of order of the clusters affect their lifetime and nucleation probability. Therefore, a single dimensional description of nucleation with a nucleus size reaction coordinate is not enough to capture the structural changes in the cluster. In addition to the crystalline nucleus size reaction coordinate, we introduce a second nucleus size variable to monitor the formation of the amorphous phase as well as the crystalline phase and show the preferred pathway of nucleation through free energy calculations. Furthermore, using a model based on previous literature (16), we study the full thermodynamic picture based on the 2D free energy surface from the molecular simulations. We then calculate the nucleation rate from the proposed mesoscopic system.

Setup of molecular simulations
We studied nucleation of NaCl from aqueous solution at two concentrations: 15.0 and 18.0 mol/kg. The NaCl interactions were calculated using the Joung-Cheatham (JC) force field (33) and the SCP/E water model (34) was used for the water molecules. The predicted NaCl solubility in water for the chosen force field is 3.7 mol/kg (30), which corresponds to a supersaturation of 4.05 and 4.8 at the chosen concentrations. The free energy calculations were done with 500 NaCl molecules and 1,851 water molecules.
All molecular dynamics (MD) simulations were done with LAMMPS Molecular Dynamics Simulator (35). Periodic boundary conditions were implemented in all directions. The simulations were run at 298 K and 1 atm. Temperature and pressure couplings were done using the Nose-Hoover thermostat (36) with a time constant of 100 time steps and Nose-Hoover barostat with a time constant of 1,000 time steps, respectively. A cut-off value of 9.0Å; was used for the Lennard-Jones and Coulombic interac-tions. Long range Coulombic interactions were calculated with the particle-particle-particle-mesh (pppm) solver with an accuracy of 10 −4 .

Selection of reaction coordinates
To monitor the formation of clusters and their structural changes, we use two collective variables: the number of ions in the largest dense cluster (n ρ ) and the number of ions the largest crystalline cluster (n c ).
The strategy for calculating the cluster sizes was inherited from Jiang et al. (25): The dense cluster size was calculated by first calculating the local density for each ion in the system. The local density was defined as the number of neighbors (ρ) an ion has within a cut-off radius r cut of 0.45 nm. An ion i was identified as "solid-like" if ρ (i) > 8. The choice of 8 for the number density cutoff was chosen based on previous literature (25), but the number density distribution of ions in a system with a crystalline slab embedded in aqueous solution supports the choice of this cut-off (see Figure S3, Supplementary Material). A clustering algorithm was used to cluster the solid-like ions, assuming that two solidlike ions belong to the same cluster if they are within 0.45 nm of each other. The freud library (37) was used to evaluate the number of neighbors and for clustering. The number of ions in the largest cluster identified with the described method was taken as n ρ . The degree of crystallinity of the ions was not considered in calculating this variable; thus, it does not carry any information about the structure of the cluster.
For calculating the crystalline cluster, the ions were first evaluated for their crystallinity by calculating their Steinhardt bondorientational order parameter (38), denoted as q 8 . The order parameter for the ith ion was calculated as follows (23,25): where q 8m (i) is given as Y 8m is the mth component of the 8th order spherical harmonics function between ion i and its neighbor j. The spherical harmonics functions were averaged over the N B = 12 nearest neighbors of the ith ion. An ion i was considered to be crystalline if q 8 (i) is greater than a threshold value of 0.45. Again, the crystalline ions were clustered assuming that two crystalline ions are connected if they are within 0.35 nm of each other. The number of ions in the largest crystalline cluster identified with the described method was taken as n c . It is worth noting that, since n ρ is calculated based only on the local density of the particles, it includes all dense particles that are crystalline as well as amorphous. Thus, for a purely crystalline particle, n c will be equal to n ρ . For a composite cluster, which is composed of crystalline particles surrounded by an amorphous layer, n c will be less than n ρ as shown in Fig. 1(a)-(c). Finally, n c can never exceed n ρ since the crystalline particles will always be dense. The choice of cut-off values and neighbor distance criteria when calculating the order parameters can result in over/underestimation of cluster sizes. This will inevitably affect our quantitative results, especially the values of the thermodynamic parameters that are obtained from fitting the composite

Free energy calculation
The 2D free energy surface as a function of the dense and crystalline nucleus sizes, βF(n), where n = (n ρ , n c ), was obtained by running 2D umbrella sampling (39) simulations with hybrid Monte Carlo/Molecular Dynamics (HMC/MD) (40). Details of the HMC/MD method can be found in the Supplementary Material along with the validation of the method in Figure S1 (Supplementary Material). The free energy is computed from where P(n) is the joint probability of observing a system with largest nucleus sizes of n c and n ρ . The value of the constant C is chosen so that βF(0, 0) is 0. The joint probability is estimated from 2D umbrella sampling simulations, using a biasing potential where K is a 2 × 2 diagonal matrix containing the biasing constants in each dimension and n * is the vector of target values to be sampled in a given simulation for n c and n ρ . It should be noted that the Landau free energy in Eq. 3 based on the largest cluster size in the system does not agree with CNT at small cluster sizes, since CNT predicts the reversible work of formation of a cluster of size n, and hence, is related to the population of nuclei of size n. The two definitions are approximately equal to each other at large clusters since it is unlikely to find two large clusters in the system (41). See the text and Figure S2 (Supplementary Material) for more detail about free energy calculations. The biasing constants were chosen as 0.15 kcal/mol and 0.07 kcal/mol for n c and n ρ , respectively, which allowed the average acceptance ratio of the MC moves to stay near 20%. The target values were varied between 0 and 50 for n c and 0 and 100 for n ρ , with enough windows to ensure that there is good overlap between samples from each window. A total of 80,000 HMC steps were collected from each umbrella window and the first 30,000 were discarded to ensure equilibrium sampling. The joint probabilities obtained from each umbrella window were combined using WHAM (42).

Mesoscopic model
The Fokker-Planck equation in multidimensional collective variable space is written as where μ is the drift term, D is the diffusivity tensor, and ρ(n, t) is the probability density with n representing phase space variables. Under equilibrium conditions, the probability flux is zero, therefore , μ(n)ρ eq (n) = ∇ · (Dρ eq (n)), In going from Eqs. 6 to 7, we use the fact that ρ eq = e −βF(n) , where F(n) is the free energy as a function of the nucleus size variables and β is 1/kT, the inverse of thermal energy. Assuming constant diffusion coefficient and substituting Eq. 7 into Eq. 5, we obtain The rate of crossing between two states can be obtained from the steady-state nonequilibrium flux of trajectories over the barrier (43,44). If the barrier is large enough, it can be assumed that the reactant side of the barrier will quickly reach equilibrium. If each trajectory that escapes to the product side is removed and replaced on the reactant side, the system will reach a nonequilibrium steady state. The steady state probability distribution of the system should satisfy the steady state Fokker-Planck equation: where ≡ ρ(n)/ρ eq (n), known as Kramer's crossover function (45).
To compute the nucleation rate constant, I, one should solve Eq. 9 and evaluate the total probability flux through a dividing surface between the solution state and the state with the nucleated cluster: where J = −ρ eq (n)D · ∇ is the probability flux, S is a dividing hyper surface in the collective variable space, and ν is the normal to that surface. The diffusivity tensor D (D ij as index notation) is evaluated by estimating the mean squared deviation (MSD) of the collective variables, n = [n ρ , n c ], from short MD trajectories and using the Einstein relation (46-48) and where δn i (t) = n i (t) − n i (t) and t f is the observed length of the trajectories. It should be noted that systematic drift was subtracted from the trajectories before they are used in the calculation of MSDs. The systematic drift, n i (t) , was evaluated by starting four independent trajectories from the same initial configuration with initialized velocities and averaging over the trajectories. This is done to prevent the contribution of any systematic drift to diffusion coefficient calculations (49). The deviations are averaged over both time and independent trajectories for improved statistics.

Results and Discussion
NaCl nucleation: minimum free energy path  Fig. 1(e) is the minimum free energy path (MFEP) on the free energy surface, which is the path that the system is free energetically most favored to follow (50). It is apparent that the dense cluster size is larger than the crystalline nucleus size at all points along the MFEP, indicating the presence of an amorphous layer around the crystalline cluster throughout the nucleation pathway. Figure 1(f) shows the change in free energy along the MFEP, in comparison with the free energy along the n c = n ρ line as well as the free energy profile obtained from the 1D umbrella sampling simulations. The MFEP includes the saddle point, which is located at n c = 22 and n ρ = 35. The comparison shows that the classical pathway requires the system to cross a higher free energy barrier compared to the MFEP, which indicates that the nonclassical pathway is thermodynamically more favored than the classical pathway at the concentration of 15.0 mol/kg. It should be noted that the MFEP still includes a single barrier, as opposed to what is expected in a two-step nucleation scenario (51,52). Therefore, we conclude that the mechanism of nucleation at this level of supersaturation is not a two-step mechanism but not necessarily classical as well. Rather, nucleation takes place with a nonclassical mechanism through the formation of composite clusters made up of crystalline and amorphous particles, namely, a composite cluster pathway. The superposition of 20 cluster configurations collected near the saddle point is shown in Fig. 2(a), which shows that the amorphous particles tend to form a layer around the crystalline cluster at the center of the cluster, in accordance with the composite cluster model.

The composite cluster pathway
The free energy change of forming a composite cluster has been formulated for nucleation of a solid from the gas phase in the presence of an intermediate liquid phase using capillarity approximation (16). The vapor-liquid-solid nucleation can be considered analogous to nucleation of a crystalline cluster from solution in the presence of an intermediate amorphous phase. To be able to compare with the simulation results, the model is reformulated here to obtain a free energy expression in terms of number of particles in the nucleus, instead of cluster radii. A schematic of the composite cluster is given in Fig. 2(b).
The three phases considered in the system are the solution phase (S), amorphous phase (A), and the crystalline phase (C). The free energy change of formation of a composite cluster in solution is the sum of the individual contributions: where F AS is the free energy of formation of the amorphous cluster in solution and F CA is the free energy of formation of a crystalline cluster within the amorphous phase. Based on the capillarity approximation, F AS and F CA can be expressed in terms of cluster sizes as where μ AS = μ S − μ A , the difference in the chemical potentials of the solution and amorphous phases. R and r are the radii of the composite and crystalline clusters, respectively, and σ ij is the surface energy of the interface between phases i and j. F CAS is a correction term to account for the short-range interaction between the C-A and A-S interfaces (16), given as where ξ is a parameter that specifies the range of interaction, and S is the spreading parameter: Equations 14-17 are combined to obtain the free energy expression Setting n c = n ρ (and R = r), Eq. 18 reduces to the CNT expression for the free energy change of single-step nucleation of a crystalline cluster from solution: from which the critical nucleus size and free energy barriers can be obtained assuming spherical geometry as where ν m is the molecular volume, assumed to be identical for the crystalline and amorphous phases. Using the dimensionless parameters from Eq. 16: Equation 18 is nondimensionalized: where Numerical values for F * CS and N * CS were obtained from fitting Eq. 19 into the free energy profile along the n c = n ρ line, shown in Fig. 1(f). The resulting parameters from the fit are shown in a Note that this value is given per mol of ion and should be doubled to obtain the chemical potential difference per mol of NaCl molecule. Table 1. The interfacial energy was determined from Eq. 20, taking ν m as the reciprocal of the number density of a NaCl crystal simulated with the JC force field (0.024 nm 3 /ion). The free energy surface obtained from simulations was scaled using the estimated F * CS and N * CS values and the nondimensional free energy expression in Eq. 22 was fitted to the scaled free energy surface. The nondimensional parameters α, β, δ, and τ are obtained from parameter estimation using nonlinear least squares. Figure 2(d) shows the resulting fit in comparison with the simulation results. It is observed that the composite cluster model agrees well with the free energy surface obtained from simulations, as the R 2 value of 0.9 suggests. The agreement between the fit and the data can be observed clearly in Fig. 2(e), where constant y lines of g(x, y) are plotted against x for both Eq. 22 and the simulation data. This shows that the nucleation mechanism for NaCl at the chosen conditions can be described well by the composite cluster model. It is worth noting that the cluster sizes involved in 0.079 ± 0.008 0.42 ± 0.05 the simulations are between 0 and 100 ions, which corresponds to cluster radii less than 1 nm 3 . At these cluster sizes, it may be expected that the capillarity approximation does not hold. The shape of the clusters is also an important factor, which is assumed to be spherical in the model. It is unrealistic to expect perfectly spherical clusters for the small sizes that this study focuses on. The sphericity analysis based on the moment of inertia tensors of the clusters showed that the clusters along the MFEP are never fully spherical but their sphericity increase as they grow, as well as their crystallinity (see Figures S6 and S7, Supplementary Material). Nevertheless, Iwamatsu's capillarity approximation-based model still agrees well with the simulation results at small cluster sizes.

Thermodynamics of nucleation pathways
Also shown in Fig. 2(d) is the transition curve. The transition curve can be defined as the collection of points located on the ridge-line between the two states, as described by Voter (53). In contrast to the 1D representation of the system where the transition state is defined by a single point, we now have a collection of nucleus sizes through which the system can cross over to the nucleated state. The transition curve gives us an idea on how the critical crystalline cluster size changes with the size of the largest dense cluster in the system. It is observed that smaller crystalline clusters are enough to initiate nucleation in the presence of larger dense clusters. The fact that the critical crystalline cluster size changes along the transition curve is in agreement with the results of Lechner et al. (54), where they employed a likelihood maximization technique to identify the order parameters that are important for homogeneous crystal nucleation in a soft core colloid model. Their results showed that information about the size of the disordered surface cloud as well as the crystalline nucleus size is necessary to predict the fate of a given nucleus. The thermodynamic parameters obtained from the fit are shown in Table 2. The negative δ value indicates that the chemical potential difference between the amorphous and the solid phase is negative. This means that the amorphous phase is less stable than the solution phase. However, σ AS , which is the surface energy of the amorphous-solution interface, is also predicted to be negative. Although unusual, negative surface free energy was recognized by Tolman (55) and reported for droplets in supercooled LJ vapor by Corti et al. (56). Surface energy dominates the thermodynamics at small cluster sizes and formation of small amorphous clusters becomes thermodynamically favorable even though the amorphous phase is less stable than the solution phase. This explains the instantaneous formation of small amorphous clusters in solution, which is observed in MD simulations of NaCl solution at 15.0 mol/kg concentration as shown in Fig. 2(c). It can be concluded that, at the concentration of 15.0 mol/kg NaCl solution, the amorphous phase is the least stable phase, followed by the solution phase and the crystalline phase, which is the most stable.

Nucleation rate
To see whether inclusion of the n ρ coordinate has an effect on the nucleation rate, we calculated the nucleation rate from the 2D system by solving Eq. 9. The diffusivity tensor that is needed for solving Eq. 9 is evaluated on the saddle point using direct MD simulations and calculating the MSD using Eqs. 11 and 12. Independent simulations were started from 25 chosen configurations that are close to the saddle point in the collective variable space. From each configuration, 4 different trajectories were started by initializing the velocities, amounting to a sum of 100 independent trajectories per point. This is done to remove any possible drift (5) as explained earlier. The resulting MSD curves and the extracted diffusivities are shown in Fig. 3(a). The components of the diffusivity tensor are estimated to be D 11 = 15.3 ± 0.9, D 22 = 3.7 ± 0.2, and D 12 = 2.6 ± 0.3 n/s. Here, the first component is the size of the dense cluster and the second component is the size of the crystalline cluster. The diffusion coefficient of the crystalline cluster is slightly less than the calculated value by Jiang et al. (25) at the same supersaturation, reported as 4.6 n/s. It is observed that the diffusion coefficient of the dense cluster is about 4 times that of the crystalline cluster, which indicates that ions can attach to the dense cluster easier compared to the crystalline cluster. This is likely because the ions need to both diffuse to the surface of the cluster and reorganize to be part of the crystalline cluster, whereas this reorganization is not needed to attach to the amorphous layer of the dense cluster.
As mentioned earlier, we assume the diffusivity to be independent of the cluster size and use the diffusion coefficients evaluated on the saddle point. This assumption is based on evaluations of diffusion coefficients at various points on the collective variable space, as shown in Fig. 3(b). It is seen that as the diffusivity of the crystalline cluster remains constant with size, there are big fluctuations in the diffusivity of the dense cluster, with no clear trend of increase or decrease with increasing size. With known free energy and diffusivity terms, Eq. 9 can now be solved to obtain the nucleation rate. Finite difference was used for numerically solving the equation with a step size of h = 0.03. Since the free energy surface is known only for discrete values of the collective variables (and the nucleus sizes themselves are discrete variables), values in between were obtained using interpolation with spline fitting. The chosen boundary conditions are shown in Fig. 3(c) and the solution for the crossover function is shown in Fig. 3(d).
The nucleation rate was obtained by evaluating the flux and using Eq. 10 as 4.3 ± 3.9 × 10 3 /s. This value is in agreement with those reported in the literature at m = 15.0 mol/kg by both brute force simulations and by using CNT rate expression (25). The agreement with CNT is remarkable, since our analysis showed that the nucleation mechanism is nonclassical at this supersaturation. There may be multiple reasons for this agreement: Although the nucleation mechanism is shown to be nonclassical in the sense that nucleation occurs through the formation of composite clusters, it is still a single-step process where the system needs to overcome a single activation barrier for nucleation. Moreover, analysis of the magnitude of the probability flux along the transition curve (Fig. 3e) shows that the majority of the probability flux occurs through a narrow distribution of nucleus sizes. This is in accordance with the classical theory that considers a single point as the transition state.

Higher supersaturation case
The nucleation mechanism and rate are analyzed up to this point at a single supersaturation of S = 4.01, which was shown to be the point at which a spinodal decomposition occurs (25). To evaluate the effect of increased supersaturation on the mechanism, we redid the free energy calculations at a concentration of m = 18.0 mol/kg, which corresponds to a supersaturation of S = 4.8. The resulting free energy surface is shown in Fig. 4(a) along with the MFEP. In comparison to the free energy surface obtained at the lower supersaturation, the MFEP is further away from the n c = n ρ line, which indicates that there is a thicker layer of amorphous particles around the crystalline core. Going from mild to severe supersaturation, the mechanism moves away from the classical pathway.
Looking at the free energy pr ofile along the MFEP (Fig. 4b), it is apparent that at m = 18.0 mol/kg, there are two distinct free energy barriers that the system needs to cross for nucleation. This agrees with the definition of the two-step process (51,52), where two barriers must be overcome. At m = 15.0 mol/kg, however, a single barrier is visible, indicating that a shift in mechanism occurs between the two concentrations. The composite cluster-free energy model given in Eq. 22 was fit onto the free energy surface obtained at m = 18.0 mol/kg and the thermodynamic parameters were obtained from parameter estimation, which are listed in Table 2. One of these parameters, μ AS , was estimated as 3.45 ± 0.25 kJ/mol ion . The positive value of this parameter suggests that the amorphous phase is now more stable than the solution phase, as opposed to the lower supersaturation case, where the amorphous phase was the least stable. The change in the sign of μ AS further confirms the shift in mechanism going from m = 15 to m = 18 mol/kg. To investigate the change in the mechanism further, we used the evaluated thermodynamic parameters in Eq. 22 and extrapolated the free energy surface to an  Table 2. The green markers denote the two saddle points on the free energy surface and the red solid lines are the MFEPs. extended collective variable space, as shown in Fig. 4(c). The extrapolation reveals the presence of a second saddle point at n ρ = 114 and n c = 1, denoted with a green marker on Fig. 4(c). At this supersaturation, there are now two possible pathways the system can follow as it nucleates. The first one is the pathway that is already visible on Fig. 4(a), where n ρ and n c increase simultaneously after the formation of amorphous clusters of about 40 ions. The second pathway, which is not visible in the free energy surface obtained from simulations, goes through the formation of large metastable amorphous clusters (n ρ > 200), after which n c starts to increase, meaning that the crystalline cluster starts to form within the amorphous cluster. The second pathway appears because the amorphous phase is now more stable compared to the solution phase, hence it is thermodynamically feasible for large amorphous clusters to form. The free energy penalties for the two MFEPs differ by 1.5 kT, with a slight thermodynamic favor for the pathway through formation of large amorphous clusters. Kinetic favor for amorphous phase nucleation was reported earlier for nucleation in undercooled liquids using density-functional simulations (57). Experimentally, Tan et al. (58) observed a decoupling between change in density and change in crystalline order in liquid to solid transition of colloids. Experimental studies of mineral nucleation report the formation of amorphous prenucleation clusters in both calcium phosphate (59) and calcium carbonate (60,61) nucleation. Using in situ TEM imaging of CaCO 3 nucleation, Nielsen et al. (60) observed multiple pathways of nucleation at varying supersaturation conditions. Our study shows that the change in mechanism could be due to a change in the stability of the amorphous clusters. Our findings also agree with the experimental observations for crystal nucleation of potassium dihydrogen phosphate (KDP) solution (62), where different pathways were observed depending on the level of supersaturation. Similar results were observed in previous computational studies of model systems. For example, a study of crystallization from solution using a Potts lattice gas model (63) showed that the preferred nucleation pathway on the cluster size-crystallinity space changed as the nucleation temperature was altered, caused by a change in relative chemical potentials of the liquid and solid phases due to the shift in temperature. The seminal paper by Ten Wolde and Frenkel (64) reported a change in nucleation mechanism from one-step to two-step as the system gets close to the liquid-liquid coexistence point for colloids with short-range attraction. Addula and Punnathanam (65) studied the two-step mechanism in crystal nucleation from Lennard-Jones vapor and reported a very similar free energy surface to Fig. 4(c), even though they examined a low driving force region of the phase diagram, at which their computed nucleation rates were virtually zero. They studied two pressures that were above the liquid-vapor coexistence curve. It could be expected that, similar to the results of this study, a further decrease in pressure into the region where the liquid phase is no longer stable would result in the disappearance of the liquid phase saddle point. A change in the free energy surface similar to the one observed in this study, depending on the relative stability of the intermediate phase compared to the mother phase was also observed for a 2D lattice model (66,67). For more realistic systems, seeing the second saddle point from all-atomic free energy calculations is challenging, since much larger box sizes would be required to avoid finite size effects. Although model extrapolation should be han-dled with caution, it gives us valuable insight about the change in the mechanism which would have been challenging to get otherwise.

Conclusions
In this paper, we examined the mechanism of NaCl nucleation from aqueous solution at elevated concentrations by providing a complete thermodynamic picture through MD simulations and free energy calculations. We showed that at S = 4.01, there is a thermodynamic preference for nucleation through the formation of clusters that have a layer of amorphous particles around the crystalline core, referred to as composite clusters. However, the free energy profile along the MFEP showed a single barrier, which indicates that nucleation occurs through a single-step mechanism, even though the mechanism is nonclassical due to the composite structure of the clusters. A free energy model that considers the formation of an intermediate phase was successfully fit against the mapped free energy profile, allowing for evaluation of thermodynamic parameters like the chemical potential differences between solution, amorphous, and crystalline phases. It was concluded that at this supersaturation, the amorphous phase is less stable than the solution phase, which explains the single-step behavior seen at S = 4.01. The nucleation rate was also computed at this supersaturation by constructing a 2D mesoscopic system. The resulting rate agrees well with the literature value obtained using CNT, which is explained by the single-step mechanism and the narrow cluster size distribution along the transition curve. As the supersaturation is further increased to S = 4.8, it was observed that the MFEP moves further away from the classical pathway. The chemical potential difference between the amorphous phase and the solution phase, obtained from the free energy model fit, revealed that the amorphous phase is now more stable than the solution phase allowing for large amorphous clusters to form in solution. The change in the relative stability of the amorphous phase resulted in the emergence of a second saddle point on the free energy surface, along with a second possible pathway for nucleation. This pathway involves formation of the amorphous clusters composed of more than 200 ions, after which the transition to a crystalline cluster takes place in a second step, in agreement with the two-step pathway.
Through a complete mapping of free energy surfaces with phase specific nucleus size coordinates at two different points on the phase diagram, we were able to demonstrate how supersaturation can change the stability ranking of the solution, amorphous and crystalline phases. The shift in the stability ranking was reflected onto the shape of the free energy surface and the preferred pathways for nucleation. Finally, we showed that the composite cluster-free energy model can effectively describe the thermodynamics not only in simple model systems (67), but also in a realistic system for solute precipitate nucleation. We believe that the methodology used in this paper can be extended to polymorphic systems to reveal thermodynamic as well as kinetic insight into competing polymorphs.

Funding
Financial support for this work was provided by AbbVie Inc. through grant numbers 8000053025 and 8000069224.