Mapping the influence of impurity interaction energy on nucleation in a lattice-gas model of solute precipitation

Abstract

We study nucleation in the two dimensional Ising lattice-gas model of solute precipitation in the presence of randomly placed static and dynamic impurities. Impurity–solute and impurity–solvent interaction energies are varied whilst keeping other interaction energies fixed. In the case of static impurities, we observe a monotonic decrease in the nucleation rate when the difference between impurity–solute and impurity–solvent interaction energies is increased. The nucleation rate saturates to a minimum value with increasing interaction energy difference when the impurity density is low. However the nucleation rate does not saturate for high impurity densities. Similar behaviour is observed with dynamic impurities both at low and high densities. We explore a broad range of both symmetric and anti-symmetric interactions with impurities and map the regime for which the impurities act as a surfactant, decreasing the surface energy of the nucleating phase. We also characterise different nucleation regimes observed at different values of interaction energy. These include additional regimes where impurities play the role of inert-spectators, bulk-stabilizers or cluster together to create heterogeneous nucleation sites for solute clusters to form.

---

1 Introduction

Nucleation is the mechanism by which a stable phase emerges from a metastable parent phase. It is the first step in the synthesis of many materials and frequently observed in nature. Classical nucleation theory^1–3 (CNT) is a well-accepted theory which can quantitatively explain this mechanism provided certain assumptions hold. A simple system of particles with short range attractive or repulsive interactions, like the 2D Ising lattice gas, exhibits nucleation behaviour accurately predicted from CNT via the Becker–Doring–Zeldovich expression.⁴ This model has been used to test assumptions and predictions of CNT in several previous studies.^5–7

In the context of solute precipitation, occupied lattice sites in this model are considered as solute particles, with empty lattice sites representing solvent. The nucleating phase transition is hence from a supersaturated solution to a precipitated solid phase. Interpreted in this way, lattice models have been used to capture two-step nucleation mechanisms⁸ and more complex nucleation pathways.⁹

Impurities or additives are often used to control the nucleation process,¹⁰ either accelerating¹¹ or decelerating^12,13 the nucleation rate. However, the mechanism by which additives influence the cluster growth varies. For example, impurities may capture the solute ions¹² or inhibit or enhance crystal growth via attachment to the surface.^11,13,14 A common real-world example where this surfactant property is used as a cleansing mechanism is the formation of micelles in aqueous solutions of soap. In recent studies, the role of structure and size of impurity particles in determining the surface properties as well as the nucleation rate is explored using molecular dynamics simulations.^15,16 It has been shown that particle size asymmetry inhibit the nucleation in spite of strong affinity between solute and impurity. In the case of polymeric impurities, the interfacial tension is decreased as the impurities bind at the interface.¹⁷ The chain length of the polymer also plays significant role in nucleation. The nucleation free energy barrier is decreased with increasing the chain length of the polymer.¹⁸ Presence of ionic impurities in aqueous nano-droplet lead to demixing of low density amorphous ice and impurity rich aqueous glass when cooled.¹⁹

Although, CNT was constructed for homogeneous nucleation in pure systems, we recently demonstrated that nucleation rates for a system containing a low concentration of homogeneously distributed impurities can be predicted by simple modification of the interfacial energy²⁰ term. In ref. 20 a special case, where both solute and solvent particles exhibit energy neutral interactions with impurities, was considered. In this scenario the impurities act as a surfactant, lowering the interfacial free energy between nuclei and the parent phase.

Studies of simple lattice models are useful to gain insight into complex nucleation mechanisms in many scenarios. Previous studies have examined the effects of surface pore-width,²¹ pore-geometry²² and surface roughness²³ on nucleation rate in the Ising model. Intelligent choice of surface geometry certainly helps to optimize the synthesis of new materials from the solution of its constituent particles. Studying competing nucleation between stable vs. metastable precipitated phases in a dimer lattice gas²⁴ provides insight into multi-component nucleation mechanism. Furthermore, examining the role of defects on magnetic droplet nucleation,²⁵ studying nucleation in the random field Ising model²⁶ and Potts lattice-gas model^18,27 are other examples where simple lattice models are have captured emergent nucleation behaviour. We note that in the Potts lattice gas model,²⁷ a surfactant-like additive was used to enhance the nucleation rate by decreasing the surface tension, similar to our previous work.²⁰

In this paper we explore a wide range of alternative scenarios, in which solute and solvent particles exhibit equal (symmetric) or opposite (anti-symmetric) interactions with impurities, over a wide range of interaction energies. We demonstrate that the surfactant-like behaviour reported previously is just one of several mechanisms by which impurities modify nucleation that can be captured within these minimal models.

The remainder of the paper is organised as follows. In Section 2, we describe the model and the algorithm used for the simulations. The impact of impurity interactions on free energy barriers to nucleation is shown in Sections 3 and 4 for static and dynamic impurities respectively. We characterize the different nucleation regimes that emerge from interaction with mobile impurities in Section 4. Section 5 contains results on nucleation rate obtained directly from forward flux sampling, which are compared to predictions made by CNT. A particular regime which exhibits clustering of impurities and enabling cross-nucleation of the stable phase is described in Section 6. Finally, we conclude in Section 7.

2 Model & algorithm

We consider a two dimensional Ising gas in the presence of randomly positioned impurities on an L × L square lattice where L = 100. Each site has a variable s_i denoting the occupancy of site i which can be either solute, solvent or impurity particles which symbolised by s_i = u, v or i respectively. We work in the semi-grand ensemble, expressing the corresponding potential aswhere ε_si,sj is the interaction energy between the species s_i and s_j when i and j are nearest neighbours, and Δμ = μ_u − μ_v is the chemical potential difference between the solute reservoir and the solvent reservoir with which the system can exchange particles. In matrix form the interaction energies can be written as:where the indices are s_i,s_j = u,v,i, starting from top left corner of the matrix. The upper-left 2 × 2 sub-matrix is the familiar Ising lattice gas with nearest neighbour interaction energy J. The remaining terms represent coupling between solute/solvent and impurity particles. We simulate below the critical temperature, such that any positive value of Δμ correspond to conditions in which the solution is supersaturated and metastable with respect to the formation of a solute rich phase.

The interaction energies of solute and solvent with impurity are denoted by ε_u,i = ε₊ and ε_v,i = ε₋ respectively throughout the paper. We study the nucleation behaviour for different values of these interaction energies at fixed impurity density. These can be categorised into three groups which are symmetric (ε₊ = ε₋), anti-symmetric (ε₊ = −ε₋) and asymmetric (|ε₊| ≠ |ε₋|) interaction energy. We set the strength of the interaction energy coupling J = 1, both within and between solute and solvent, and the impurity–impurity interaction strength is set to 0 throughout the models studied in the paper.

We simulate this model using Metropolis Monte Carlo dynamics. The usual transmutation moves of solute into solvent and vice versa represent removal of one species to its corresponding reservoir and replacement from the other reservoir.

In addition we introduce moves which allow the impurities to migrate. We define a mobility parameter α of the impurities such that α = 0 corresponds to static impurities and α = 1 corresponds to the fastest moving impurities with no transmutation dynamics. At each Monte Carlo move we randomly generate a random number ξ uniformly distributed between 0 and 1. If ξ < α, we attempt a “non-local swap” move in which the occupancy of a randomly selected impurity site is swapped with a second randomly selected site a distance d away. Specifically we generate a displacement vector (Δx,Δy) with circular symmetry by setting Δx to a random integer between −d to d then define . Otherwise if ξ ≥ α we attempt to transmute a randomly selected solute or solvent. One Monte Carlo step consists of L² moves of any type and each move is accepted according to the usual Metropolis acceptance criterion.

We set the linear distance between two sites involved in a swap move to d = 4 with the nearest integer approximation. The non-local swap move of impurities is introduced for fast equilibration of the system. The impurities can leave or enter the cluster efficiently with implementation of these moves. We note that this non-local impurity swap dynamics is different to the local (nearest neighbour) impurity swap dynamics studied in ref. 20, and so absolute estimated rates for the same α should not be compared directly.

For constant d the parameter α controls the mean squared displacement of the mobile impurities with time, which increases monotonically with α. We do not tune this parameter to match any particular real system, rather we explore the two limiting cases of static and “fast” impurities. The latter implies the spatial distribution of impurities equilibrates on timescales which are rapid compared to changes of cluster size. This could represent any system in which cluster growth is attachment limited and mobile impurities are present.

We analyze the nucleation behaviour of the system by studying the nucleation free energy barrier for different impurity interaction energies. Since the formation of a post-critical cluster from the metastable phase at low temperatures is a rare event, unbiased simulations fail to sample the configuration space sufficiently within a tractable computational timescale. Therefore we use the umbrella sampling^28,29 method to compute the nucleation free energy barrier as a function of solute cluster size. We use the standard geometric definition of a cluster, i.e. a set of solute particles that are contiguously connected via nearest neighbours to other solute particles. An infinite square-well potential spanning a cluster size range of 20 is used to bias the system to remain in a particular cluster size window, and we simulate multiple overlapping windows to cover the full range of relevant cluster sizes. The segments of the free energy curve obtained within each window are then combined with the appropriate shift to reconstruct a smooth and continuous free energy curve. For windows which span smaller cluster sizes it is possible for multiple clusters to appear simultaneously that satisfy the size criterion of that window. It is important to count all such clusters and not just the largest in order to construct a free energy barrier consistent with CNT and the Becker–Doring–Zeldovich nucleation rate calculation. Otherwise the free energy exhibits a spurious minimum at the most frequently occurring largest cluster size observed in the metastable parent phase, which is grid size dependent.

The nucleation rate, i.e. the rate per unit area of forming post-critical nuclei from the metastable solution, is another important quantity that can be calculated either by direct simulation, or by using classical nucleation theory. Forward flux sampling^30–33 is a direct simulation method used for calculating the nucleation rate. In the forward flux sampling method, we define a set of “interfaces” at increasing values of the largest cluster size λ. The nucleation rate from the metastable solution phase can be written as

where I₀ is the positive flux (crossings per unit time) measured through the zeroth interface in an unbiased simulation of the metastable phase. The quantity P(λ_i+1|λ_i) is the probability of a simulation initialised at the i-th interface reaching the (i + 1)-th interface before returning to the metastable solution. The interface λ₀ is chosen to be some small cluster size such that sufficient sampling of crossings can be sampled within a tractable simulation time. Use of the largest cluster size does result in a value of I₀ which is system size dependent for a given choice of λ₀, however the overall rate is not.³⁴

The details of our implementation of both umbrella sampling and forward flux sampling algorithms could be found in ref. 20.

3 Static impurities

We wish to explore the full range of behaviour as a function of ε₊ and ε₋. In the case of static impurities only transmutation moves are performed. The interaction energies ε₊ and ε₋ appear in the acceptance probability for these moves via their difference ε_d = ε₊ − ε₋ since the moves replace impurity-solute interactions with an equal number of impurity-solvent interactions or vice versa. In Fig. 1 we have pictorially demonstrated the energy required in transmuting a solvent to a solute with one and two surrounded impurities. For static impurities is hence possible to reduce the exploration of interaction energies to a single parameter ε_d. For convenience of comparison with later results, we choose to set ε₊ = +ε and ε₋ = −ε such that ε_d = 2ε. This sets impurity–solvent and impurity–solute interactions to be anti-symmetric in our simulations, however results for static impurities would be numerically identical for any choice of ε₊ and ε₋ that preserves ε_d.

Fig. 1 The energy required to transmute a solvent (white) to a solute (red) surrounded by one and two impurities (blue) for the configurations depicted in (a) and (b) are respectively δε = ε_d − 2J and δε = 2ε_d − 4J. It depends only on ε_d = ε₊ − ε₋ and J. The configurations before and after updating the central site are shown in left and right side of the arrow.

We calculate the nucleation free energy βF(λ) as a function of cluster size λ over a range of the dimensionless impurity interaction strength βε and at several values of static impurity density. Here β = 1/k_BT is the inverse temperature T, and k_B is the Boltzmann constant which we set to 1. In the free energy calculation, configurations in each umbrella sampling window are sampled over 48 realisations of the static impurity disorder and combined to compute a single free energy profile. We then numerically estimate the critical cluster size λ_c and free energy barrier height βF(λ_c) for different values of βε_d/2. The free energy obtained from the umbrella sampling calculations is fitted by the free energy expression^4,20

where B = −log ρ₁ − A_s + A_b, and ρ₁ is the density of isolated solute particles in the solution phase. The bulk term A_b = βΔg, where Δg is the bulk free energy difference per particle between solute particles in the stable nucleating phase and in the metastable solution phase. At low temperatures, where the two phases are dominated by solute and solvent respectively, Δg ≈ Δμ. Estimation of λ_c and βF(λ_c) is done by fitting βF(λ) to the expression given in eqn (4) using A_b and A_s as fitting parameters and finding the position and value of its maxima. However, if the free energy is not well fitted by eqn (4), we add a higher order polynomial term, i.e., λ^3/2 with its prefactor as fitting parameter to improve estimation of λ_c. The higher order term is required for static impurities when βε_d is large and positive. Here the free energy shows significant deviation from the expression shown in eqn (4) due to the strong repulsion and confinement effect from the impurities.

The values of βJ we study are chosen to be greater than the critical inverse temperature of the two dimensional Ising model .³⁵ However, at high values of βJ, when the corresponding temperature is low, the rare-event sampling algorithms used in this paper become inefficient. We hence focus on intermediate values of βJ which are 0.67 and 0.83, while doing umbrella sampling and forward flux sampling simulations. We note that performing simulations at different values of βJ implies performing it at different temperatures as the coupling constant is fixed to J = 1. We choose different values of βJ in simulations to show the validity of our results for a range of temperature.

Variation of the free energy barrier height and the critical cluster size as a function of interaction energy difference βε_d, for ρ_i = 0.004 (green squares) and 0.008 (blue up-pointing triangles) at βJ = 0.67 and βΔμ = 0.067, with static impurities are shown in Fig. 2(a) and (b) respectively. The barrier height as well as critical cluster size increase with increasing βε_d.

Fig. 2 Free energy barrier height (a) F(λ_c) and (b) critical cluster size λ_c as a function of the dimensionless magnitude of impurity interaction energy βε for βJ = 0.67 and βΔμ = 0.067. Static (α = 0) and dynamic (α = 0.05) impurity cases with densities ρ_i = 0.004 and 0.008 are plotted. In both cases anti-symmetric impurity interactions (ε₋ = −ε₊) are used. The free energy barrier height tends toward saturation with increasing βε for ρ_i = 0.004 in the case of both static (α = 0) and dynamic (α = 0.05) impurities denoted by green squares and red circles respectively. At higher impurity density (ρ_i = 0.008), the saturation is also observed in the case of dynamic impurities as denoted by magenta down-pointing triangles, but not in the case of static impurities as denoted by blue up-pointing triangles. Similar behaviour is also observed for λ_c as a function of βε. Dashed horizontal lines represent the barrier height βF(λ_c) = 40.81 and critical cluster size λ_c = 496 for the pure system (ρ_i = 0) at βJ = 0.67 and βΔμ = 0.067.

This trend in free energy barrier is expected since the energy required to transmute a solvent to a solute (in the presence of a neighbouring impurity) is lowered when βε_d is negative, meaning that nucleation can proceed preferentially in regions where impurities are present. For large negative βε_d, we would expect to see spinodal decomposition as a result of strong attraction between the impurities and solute.

For positive βε_d nucleation will preferentially occur away from impurities and so the nucleation rate will saturate with respect to βε_d provided there is sufficient space for a critical nucleus to form without needing to neighbour any impurity sites. This seems to be the case for ρ_i = 0.004, but not for the higher impurity density of ρ_i = 0.008 where the barrier height continues to increase with βε_d indicating that critical nuclei cannot form without encountering impurities that impede their formation.

In Fig. 2(b), critical cluster size vs. βε_d plots show similar behaviour, however growth is faster compared to the free energy barrier. See Fig. S3 in the ESI† for detailed free energy plots.

The barrier heights and critical sizes are consistent with our previous study when ε_d = 0.²⁰

We note that our free energy curves are computed by sampling configurations in each λ window over several realisations of static disorder. Calculations at large and positive βε may be dominated by a small number of these realisations where sufficient room is available to form a critical nucleus without encountering impurities. Such clusters would have low energy and hence higher probability of formation when compared to equal size clusters in other realisations of the impurity disorder. This observation is related to that made in a recent study of the 3D random field Ising model²⁶ in which a spatially dependent reaction coordinate is used to account for the position, as well as the size of a nuclei, since the random field can create preferential locations for nuclei to form.

4 Dynamic impurities

In the case of dynamic impurities we include non-local impurity swap moves as described in Section 2. The change in potential energy δΦ resulting from these moves depends on J, ε₊ and ε₋ explicitly and cannot be reduced to fewer parameters as in the static case. This is because impurities can hop between two sites that have differing numbers of solute and solvent neighbours. We study the system for different interaction types, e.g., symmetric (ε₋ = ε₊), anti-symmetric (ε₋ = −ε₊) and asymmetric (|ε₋| ≠ |ε₊|) with respect to impurity interactions with solute and solvent particles. We define the mobility parameter α, which can take values between 0 ≤ α < 1. The case α = 0 represents static impurities and α ≠ 0 represents dynamic impurities with mobility increasing with α. For the non-local impurity swap moves used in this study, equilibration of the impurity distribution is fast compared to the timescale on which clusters of solute particles grow or shrink. The free energy curves presented below hence represent a thermodynamically controlled nucleation process in which the impurity the distribution is sampled from a quasi-equilibrium distribution at each cluster size.

4.1 Anti-symmetric interaction energy

Here we set the interaction energies ε₊ = −ε₋ = ε. Positive values of ε make impurity–solute interactions unfavourable and impurity–solvent interactions favourable. For negative values of ε the preference is reversed. A sequence of free energy curves as a function of cluster size λ for different interaction energies βε with βJ = 0.83, βΔμ = 0.083, ρ_i = 0.02 and α = 0.05 is shown in Fig. 3. The impurities act as nucleating sites for βε < 0 (lowering free energy barrier height), and repel solute for βε > 0 (higher free energy barrier height).

Fig. 3 Nucleation free energy with dynamic impurities, varying anti-symmetric interaction energy βε₊ = −βε₋ = βε for βJ = 0.83, βΔμ = 0.083, ρ_i = 0.02 and α = 0.05. Dotted lines are obtained from numerical fit of eqn (4) varying bulk A_b and surface A_s terms.

The trend in barrier height and critical cluster size with βε is compared to the static impurity case in Fig. 2. For large positive values of βε we expect the formation of clusters to exclude impurities even at high impurity density due to migration of impurities away from growing clusters (not possible for the static impurities). At impurity density ρ_i = 0.004 and βε ≳ 0.53, the impurity interaction energy with solute particles is sufficiently unfavourable that all impurities are excluded from the growing cluster and there is no further change in barrier height as shown in Fig. 2(a) by red circles. The saturation threshold is βε ≳ 0.8 when ρ_i is increased to 0.008, denoted by magenta down-pointing triangles. Similar saturation is seen in critical cluster size as shown in Fig. 2(b). See Fig. S1(a) in the ESI† for complete free energy plots.

A similar saturation in free energy barrier is observed even at impurity density ρ_i = 0.02 and a different βJ = 0.83, but now at a more positive βε ≳ 1.17. See Fig. S1(b) in the ESI† for free energy plots.

Unlike the static impurity case, this saturation of the barrier height with increasingly positive βε results in a lower bound on the nucleation rate when the impurities are dynamic as the nucleation barrier height cannot increase further. As with the static case, there is no equivalent upper bound on the rate with dynamic impurities, as decreasingly negative βε will lower the free energy of any clusters containing impurities until the spinodal limit is reached and the system can spontaneously transform into the solute-rich phase.

As the interaction strength βε is made increasingly positive, we see a crossing both in βF(λ_c) and λ_c when these quantities are plotted for two different impurity densities, ρ_i = 0.004 and 0.008. This occurs at the value βε* ≈ 0.13 for βJ = 0.67, βΔμ = 0.067 and α = 0.05 (see Fig. 2). Notably this is the value of βε for which the system exhibits near identical barrier height and critical cluster size to the case where no impurities are present, implying cancellation of competing effects.

The role of impurities is opposite either side of this crossover. In region βε > βε*, when impurities repel solute particles weakly, the barrier height is increased with increasing ρ_i. This is analogous to an observation in ref. 13, where the growth of succinic acid is inhibited by different additives (glutaric acid, heptanedioic acid and azelaic acid), effectively increasing the interfacial tension with increasing additive concentration. For βε < βε*, when there is a weak interaction between impurities and solute, interfacial tension decreases with increasing ρ_i enhancing the nucleation rate. This is analogous to the experimental observation in ref. 11, where the presence of type-III antifreeze protein enhance the growth of ice nucleation by sitting at the boundaries of the cluster. A similar behaviour is obtained in a simulation of Potts lattice-gas in the presence of low dosage additives.²⁷ The two curves cross where these competing effects cancel.

In Fig. 2, while comparing the dynamic (α = 0.05) vs. static (α = 0) plots, the barrier height and critical cluster size for dynamic impurities is always less than the values corresponding to static impurities, independent of interaction energy βε. This implies that the faster nucleation rate compared in the static case does not depend on the repulsive or attractive nature of the microscopic interactions when impurities are mobile. This behaviour is also observed for neutral impurity interactions in ref. 20.

4.2 Symmetric interaction energy

Now we set the interaction energy ε₊ = ε₋ = ε to be symmetric. In this case the total impurity-solute and impurity-solvent interaction energy takes the constant value 4ερ_iL², independent of the impurity location, provided the impurities are at sufficiently low density that impurity–impurity interactions are negligible, which we expect to be the case for negative βε. The quantity ρ_iL² is the total number of impurity particles present in the solution. This impurity interaction term in the total energy can hence be treated as constant energy shift in the total energy which makes no difference to the nucleation behaviour. Thus we expect the free energy as well as the nucleation rate to be unchanged with respect to varying symmetric interaction energy ε as long as the impurities are sparsely distributed. In Fig. S5 of the ESI,† we have plotted the variation of free energy for different negative values of symmetric interaction energy between impurity–solute and impurity–solvent, for βJ = 0.83, βΔμ = 0.083 and α = 0.05. We do not see a significant variation in free energy barrier for the range of (negative) βε plotted in Fig. S5 (ESI†) suggesting no impurity clustering occurs over that range. For positive βε however, clustering of impurities will be expected, reducing the number of unfavourable impurity–solvent and impurity–solute interactions. This will be explored more generally in the next section.

4.3 Behaviour map for asymmetric interactions

If allowing the impurity–solute and impurity–solvent interaction energies to be neither symmetric nor anti-symmetric, the available nucleation behaviour is richer. We aim to map the possibilities by characterising the nucleation behaviour at each interaction choice. It is expected that for different interaction energies ε₊ and ε₋ the impurities preferentially occupy different positions in relation to the growing solute nucleus. These include impurities completely inside the clusters, completely outside the clusters and at the boundary of the solute clusters. To construct a map of this behaviour we calculate the average fraction of impurities that are located at the boundary of the largest cluster at fixed βε₊, βε₋ and βΔμ. The local micro-state of impurities can be divided into five different groups depending on the number of nearest neighbour solutes (0 to 4). We count the fraction of impurities ϕ that have one, two or three nearest-neighbour solute particles as shown in Fig. 4. This count takes place within a biased simulation in which we restrict the size of the largest cluster to be between 800 to 1000, and the size of the next largest cluster to be less than 30% of the largest cluster size to avoid contacting with the largest cluster when the nucleation rate is high. We count only the impurity sites attached to the largest cluster and plot ϕ as a function of βε₊ and βε₋. The resulting map for βJ = 1.11, βΔμ = 0.11, ρ_i = 0.02 and α = 0.05 is shown in Fig. 5(a). We have used a slightly higher value of βJ here as the preferential occupancy of impurities at the boundary of the growing cluster is only observed at low temperatures or equivalently at high βJ.²⁰ With decreasing βJ the intensity of bright area decreases, making it difficult to differentiate between four regimes of the behaviour map introduced in next paragraph.

Fig. 4 Three types of impurity micro-states that contribute to the accumulation of impurities at the boundary of the nucleus. Red and blue boxes represent impurity and solute respectively. A shaded box could be either solvent or impurity.

Fig. 5 (a) Variation of average impurity density at the boundary of the largest cluster ϕ as a function of βε₊ and βε₋ with dynamic impurities at density ρ_i = 0.02, βJ = 1.11, βΔμ = 0.11 and α = 0.05. Depending on the positional occupancy of the impurities the interaction energy space could be divided into four regimes (1) surfactant: impurities prefer to occupy the boundary of the cluster, (2) inert spectator: impurities are completely excluded from the cluster without taking part in the nucleation process, (3) insoluble heterogeneous nucleation sites: impurities form clusters which act as nucleation sites and (4) bulk stabilizer: impurities are completely inside the clusters stabilizing the bulk phase. Snapshots of the system with biased simulations at quasi-equilibrium with largest cluster size λ bounded between [800, 1000] at regime (b) 1, (c) 2, (d) 3 and (e) 4 of the behaviour map.

The behaviour map can be divided into four regimes depending on the positional occupancy and role of impurities. These regimes are (1) surfactant: the bright area where impurities prefer to occupy the boundary positions of a cluster acting as a surfactant, (2) inert spectator: the blue area in the right side of the behaviour map, where impurities are excluded from the nucleating clusters of solute, (3) insoluble heterogeneous nucleating sites: the noisy bright area at the top-right corner, where the impurities themselves form clusters which can act as heterogeneous nucleation sites (see Section 6), and finally (4) bulk stabilizer: the blue area in the left side of the behaviour map, where impurities are preferentially located inside the red clusters as inclusions, stabilizing the bulk phase. Approximate boundary lines between regimes are drawn by black dotted lines in Fig. 5(a). The previously discussed symmetric and anti-symmetric cases correspond to behaviour along the two diagonals. Snapshots of configurations from each of the four regimes in quasi-equilibrium with largest cluster size λ confined between [800–1000], are shown in Fig. 5(b)–(e) for different interaction energies (b) βε₊ = −1.67, βε₋ = −1.67 (surfactant), (c) βε₊ = 1.67, βε₋ = −1.67 (inert-spectator), (d) βε₊ = 1.67, βε₋ = 1.67 (insoluble impurity clusters) and (e) βε₊ = −1.67, βε₋ = 1.67 (bulk-stabilizer).

An ad hoc way of estimating the extent of regime (1) (impurities as surfactants) is the following. Impurities prefer to be in the solution phase and solute phase when ε₊ > ε₋ and ε₊ < ε₋ respectively. Arguably, ε₊ = ε₋ would correspond to the regime where impurities prefer to occupy the boundary of a cluster. The expected width of region (1) on either side of this diagonal would be obtained by analyzing the stability of the interface. If we consider a flat interface without impurities separating pure solute and pure solvent regions, the interface energy per unit length would be J. On the other hand, if we add one layer of impurities at the interface between solute and solvent phase, the interface energy per unit length would be (ε₊ + ε₋)/2. The interface with impurities would be stable if the condition (ε₊ + ε₋)/2 < J is satisfied. Combining these two criterion, we obtain βε < βJ, where ε₊ = ε₋ = ε. We see that the derived condition for the width of surfactant regime is approximately satisfied in Fig. 5(a).

We analyze the impact of impurities in the regime outside the surfactant area and relate these to the trends in free energy observed in earlier sections. In the surfactant regime, both solute and solvent attract impurities and we see surface accumulation of impurities. However, the free energy barrier does not change significantly with varying symmetric interaction energy. In the inert-spectator regime, solute repels but solvent attracts the impurities and we see strong exclusion of impurities from the nucleus. In this regime the barrier height remains unchanged with increasing symmetric interaction since all impurities are excluded from the nucleus. The regime in which impurities act as heterogeneous nucleation sites, both solute and solvent are repelled by the impurities. This strong repulsion forces impurities to form clusters. In the bulk-stabilizer regime, solute attracts but solvent repels impurities making nucleation strongly favourable with the impurities acting as nucleants. Here we see the presence of multiple clusters with impurities as inclusions. We also see a low barrier height to nucleation in this regime.

It is evident from Fig. 5(a), that the impurities preferentially occupy the boundary positions of the cluster when ε₊ ∼ ε₋. But there exists a small asymmetry between these two energies as the maxima of the bright regions in Fig. 5(a) does not go exactly through the diagonal, i.e., the βε₊ = βε₋ line. To illustrate this asymmetry we change the co-ordinate system of our behaviour map from (βε₊,βε₋) to (βε₊ + βε₋,βε₊ − βε₋). Projections of the transformed behaviour map along β(ε₊ + ε₋) = c line for different constant values of c are shown in Fig. 6. The maxima of the projection plot occurs at β(ε₋ − ε₊) ≈ 0.51 for βJ = 1.11 and βΔμ = 0.11. This implies impurities should have a small energy bias towards the solute compared to the solvent to maximise surface accumulation when the largest cluster is of the particular size used to construct this map. We expect the presence of a transition with decreasing βJ (equivalent to increasing temperature when J is fixed) to a situation where impurities would no longer act as a surfactant for any type of interaction energy. This was analyzed for neutral impurities in ref. 20.

Fig. 6 One dimensional projection of the average impurity density at the boundary of the largest cluster, ϕ along β(ε₋ + ε₊) = c line for different c at fixed impurities of density ρ_i = 0.02 with βJ = 1.11 and βΔμ = 0.11. The maxima occurs at β(ε₋ − ε₊) ≈ 0.51. The position of maxima is shown by vertical solid line.

5 Nucleation rate & Becker–Doring–Zeldovich analysis

We use forward flux sampling (FFS) to calculate the nucleation rate, i.e., the rate at which post-critical solute clusters that reach macroscopic size are obtained from the initial metastable solution phase. The mathematical expression for the nucleation rate within FFS is given in eqn (3). The right hand side of eqn (3) may be interpreted as the rate of obtaining a cluster of size λ = λ_n+1 at the (n + 1)-th interface, from the solution phase. The lower limiting value at which the decreasing rate converges for λ > λ_c is the nucleation rate. We denote L × L Monte Carlo moves as the unit of time, i.e. one Monte Carlo step. In the case of static impurities, the nucleation rate I is measured by the crossings per unit Monte Carlo step per single site which is consistent with the definition use by other authors.^4,25 For dynamic impurities, we divide the number of crossings per unit Monte Carlo step per single site by (1 − α), i.e. time is only progressed by attempted transmutation moves between solute and solvent. Impurity dynamics are considered fast on this timescale such that the time elapsed during the non-local swap moves can be neglected.

The rate I(λ) of obtaining a cluster of size λ starting from the solution phase is plotted in Fig. 7 for impurity density ρ_i = 0.012 at βJ = 0.83, βΔμ = 0.083 and α = 0.05. Curves are plotted for a range of values of βε that lie on the diagonal line joining the top-left (bulk stabilizer) to bottom-right (inert spectator) corner of the behaviour map [see Fig. 5(a) for the behaviour map with impurity density ρ_i = 0.02 at βJ = 1.11, βΔμ = 0.11 and α = 0.05]. The nucleation rate does not further decrease beyond a minimum value for βε ≳ 1.33 which belongs to the inert-spectator regime of the behaviour map. This represents the limit beyond which the probability of finding an impurity inside, or at the boundary of a large solute cluster is negligible and hence there is no further impact on the nucleation rate with increasing impurity interaction energy.

Fig. 7 Rate of obtaining a cluster of size λ from the metastable solution phase for different anti-symmetric interaction energies βε = βε₊ = −βε₋ at fixed βJ = 0.83, βΔμ = 0.083 and ρ_i = 0.012 with dynamic (α = 0.05) impurities.

The nucleation rates extracted from Fig. 7 are re-plotted in Fig. 8(a) (red curve) and compared to the rates along the opposite diagonal of the behaviour map, where impurity interaction energies are symmetric (blue curve). Here the lower limit of the nucleation rate for anti-symmetric impurity interactions is visible as saturation of the rate I with respect to increasingly positive βε. It can also be seen that the nucleation rate increases without apparent limit as βε becomes more negative. This represents tending toward spinodal decomposition as the stability of the bulk solute phase is enhanced by the presence of impurity inclusions, similarly to the static case.

Fig. 8 (a) Nucleation rates as a function of interaction energy βε for anti-symmetric and symmetric interaction energy at fixed βJ = 0.83, βΔμ = 0.083, ρ_i = 0.012 and α = 0.05. (b) Saturated (with respect to increasing βε) nucleation rate I_s(ρ_i) as a function of ρ_i with anti-symmetric interaction energy with fixed βJ = 0.83, βΔμ = 0.083 and α = 0.05. The estimated standard error in calculating the nucleation rate plotted both in (a) and (b) is less than the size of the symbols.

A natural question to ask is whether the minimum rate on increasing βε for anti-symmetric interactions depends on impurity density. To answer that, we plot the saturated nucleation rate I_s(ρ_i) for large and positive βε and anti-symmetric interactions at different impurity density ρ_i as shown in Fig. 8(b). We see that the saturated rate increases monotonically with decreasing ρ_i. As the impurities are excluded from the cluster all of them enter into the solution. The impurity density in the solution increases with increasing ρ_i. This excess impurities present in the solution change the effective chemical potential difference Δμ between the solution phase and the crystalline phase when ρ_i is varied, changing the nucleation rate. We note that in the inert spectator regime while going from ρ_i ≠ 0 (three-component solution) to ρ_i = 0 (two-component solution), the effective Δμ also changes which has a slight impact on nucleation rate.

In the case of symmetric impurity interaction energies, the two ends of the (blue) curve in Fig. 8 lie in the surfactant regime (negative βε) and in the regime where impurities act as heterogeneous nucleation sites (positive βε). As expected from our analysis in Section 4.2, we do not see a significant variation in the nucleation rate for symmetric interactions where βε is negative and impurity–impurity interactions are rare. For large and positive βε the limiting behaviour is that of a single large impurity cluster with a surface energy independent of whether it is surrounded by solute or solvent. Examination of simulation snapshots shows that the intermediate regime (βε ≈ 1) is characterised by the presence of both isolated impurities surrounded by solvent, and a single substantial impurity cluster. Nucleation of solute clusters occurs preferentially at the interface between this impurity cluster and the solution, but it is unclear why the nucleation rate is slightly enhanced in this case compared to larger values of βε where all impurities are present within a single cluster.

Becker–Doring–Zeldovich analysis: Studying the ability of a Becker–Doring–Zeldovich analysis^4,20 to reproduce trends in nucleation rate with varying impurity interactions can be instructive. In particularly it gives an understanding of which physical parameters (surface versus bulk free energies, kinetics) must be varied to fit the numerical simulation data and hence verify our mechanistic interpretation of results in the observed regimes.

We fit the free energies obtained from the umbrella sampling calculations with the modified free energy expression given in eqn (4). In the case of dynamic impurities, we use A_b and A_s as fitting parameters to fit eqn (4) with the free energy curves obtained using umbrella sampling simulations, shown in Fig. 3 for ρ_i = 0.02. The values of A_b and A_s obtained from the fitting is shown in Table 1. These are used to calculate the nucleation rate described in the next paragraph. We see monotonic increase and decrease in A_b and A_s respectively with decreasing βε.

Table 1 Estimates of A_b and A_s from fitting eqn (4) to the free energy plots obtained from umbrella sampling simulations shown in Fig. 3 for βJ = 0.83, βΔμ = 0.083, ρ_i = 0.02 and α = 0.05 with anti-symmetric interaction energy and dynamic impurities with α = 0.05

βε	A b	A s
1.33	0.056	4.54
0.83	0.057	4.47
0.17	0.069	4.16
0	0.082	3.81
−0.17	0.114	3.10

---

The Becker–Doring–Zeldovich expression of the nucleation rate can be written as

I_BDZ = D_cΓe^−βF(λ_c), (5)

where D_c is the diffusion coefficient, Γ is the Zeldovich factorevaluated at the critical cluster size λ_c. Using eqn (4), these quantities can be expressed asThe diffusion coefficient D_c = 〈Δλ(t)²/2t〉 is obtained from independent simulations, starting from the critical cluster size and computing the slope of the mean square displacement versus time. See Fig. S9 in ESI† for D_c estimation. In Table 2, we compare the nucleation rates obtained from the Becker–Doring–Zeldovich analysis (see eqn (5)) and forward flux sampling simulation for different anti-symmetric interaction energies βε at dynamic impurity density ρ_i = 0.004. The final two columns, I_BDZ and I_FFS, are the rates obtained from eqn (5) and independent forward flux sampling simulations respectively. The maximum error in determining I_BDZ is obtained by the expressionwhere Δx is the error in determining the quantity x. The excellent agreement between the results within the calculated error validates the application of CNT for the model studied, with only refitting of the surface and bulk terms due to the presence of impurities required. See Tables S1 and S2 in the ESI† comparing I_BDZ and I_FFS at different regimes of the behaviour map.

Table 2 Comparison of nucleation rates obtained from Becker–Doring–Zeldovich analysis (I_BDZ) and forward flux sampling (I_FFS) for βJ = 0.67, βΔμ = 0.067, ρ_i = 0.004 and α = 0.05 with anti-symmetric interaction energy. The maximum error in determining I_BDZ and I_FFS are 72% and 10% respectively

βε	λ c	βF(λc)	D c	I BDZ	I FFS
0.8	564	42.94	41.2	2.8 × 10^−20	2.5 × 10^−20
0.13	510	40.6	39	3.1 × 10^−19	3.1 × 10^−19
0	469	38.8	37.1	1.8 × 10^−18	2.4 × 10^−18
−0.13	388	35.2	33.7	6.7 × 10^−17	1.0 × 10^−16

---

6 Impurity clustering and cross-nucleation

Impurities form multiple clusters when both impurity–solute and impurity–solvent interaction energies are positive, as seen in [Fig. 5(d)]. Within this regime, when the repulsive interaction with solute and solvent is sufficiently strong, a single impurity cluster becomes the most stable configuration. We observe that nucleation of the solute phase starts from the boundaries of the impurity cluster, although the interaction between impurity and solute is repulsive.

Snapshots of a nucleating system, with symmetric interaction energy βε = 1.33, βJ = 0.83, βΔμ = 0.083, ρ_i = 0.012 and α = 0.05, obtained from the umbrella sampling simulation for window-15 (largest cluster size lies between 150 and 170) and window-104 (largest cluster size lies between 1040 and 1060) are shown in Fig. 9(a) and (b) respectively. The binding between impurity and solute clusters occurs because the total surface energy of the two clusters is reduced when impurity–cluster and solute–cluster share a common boundary compared to the case when they are separated.

Fig. 9 Snapshots from the umbrella sampling simulations, at (a) window-15 (solute cluster size lies between 150 and 170) and (b) window-104 (solute cluster size lies between 1040 and 1060), showing the binding of an impurity clusters and solute cluster for symmetric interaction strength βε = 1.33, βJ = 0.83, βΔμ = 0.083, ρ_i = 0.012 and α = 0.05.

Assuming a circular shape of both nuclei we may write the surface energy difference between the bonded-cluster and two separate-clusters as σ_b(R,r) − σ_s(R,r) ≈ −2rJ + πr(ε₊ − ε₋) (see Appendix A for derivation), for R ≫ r, where R and r are the radius of solute–cluster and impurity–cluster respectively. For symmetric interaction energy ε₊ = ε₋, the second term in the right hand side vanishes and the surface energy difference becomes completely negative stabilising the bonded configuration. This is an example where the attraction between two clusters is induced by the microscopic repulsion between two particle types. A similar idea has been used to calculate the free energy of a droplet starting to grow at the boundary of a parent nucleus.³⁶

In the current context, this preferential formation of solute clusters at the boundary of impurity clusters can be considered as cross-nucleation.^37,38 Here, the impurity cluster acts as a heterogeneous nucleation site for the nucleation of solute clusters.

7 Conclusion

We have studied the nucleation behaviour of a two dimensional Ising lattice-gas model in the presence of static and dynamic impurities with varying impurity–solute and impurity–solvent interaction energy.

In the case of static impurities, we have shown that the nucleation free energy barrier height increases on increasing the difference between impurity–solute and impurity–solvent interaction energy ε_d. The barrier height shows saturation with increasing ε_d when the static impurity density is low. However, we do not see such barrier height saturation when impurity density is high enough so that a critical cluster cannot fit into the largest void space between impurities.

In the case of dynamic impurities, at high βJ (or equivalently, low temperatures) we observe preferential occupancy of the impurities at the boundary positions of the nucleus when the interaction energy of impurities with solute and solvent are similar. We have studied the system with varying the interaction energy and characterised four different nucleation regimes depending on the role and positional occupancy of impurities in the nucleation process. These regimes are surfactant, inert spectator, heterogeneous nucleation sites of impurity clusters and bulk stabilizer. Free energy behaviour and nucleation rate have been studied in each regime and the limits of impurity influence have been established.

In this paper, we have the interactions between impurities to be neutral. Given the non-trivial behaviour when impurity clusters form, it would be interesting to extend this work for non-zero impurity–impurity interaction energy. How the different regimes in the behaviour map change with varying impurity–impurity interaction would also be interesting to investigate. It might be argued that the Monte Carlo moves in our current model limit the study to regimes where certain kinetic assumptions apply. It may be interesting to extend this to include diffusive transport of solute and solvent in place of transmutation moves, allow of concentration gradients of impurity, and other modifications to determine if new regimes of nucleation behaviour emerge.

Author contributions

DM developed the necessary simulation software, performed the simulations and wrote the manuscript. Both authors contributed to the design and interpretation of the simulations. DQ edited the manuscript.

Data availability

Data associated with this manuscript is available via the University of Warwick Research Archive Portal at https://wrap.warwick.ac.uk/187431/.

Conflicts of interest

The authors have no conflicts of interest to declare.

Appendices

A Stability of bonded impurity–cluster and solute–cluster

A schematic representation of a bonded impurity cluster and solute cluster is shown in Fig. 10. The surface energy of the bonded configuration σ_b(R,r) can be written as the surface energy contribution obtained from the green boundary line which can be expressed as

σ_b(R,r) ≈ (2πR − 2r)J + πr(ε₊ + ε₋), (10)

when we assume circular shape of the clusters. We also assume that radius of the solute–cluster does not change for the separate configuration which is true when R ≫ r. Now, the total surface energy for separate configurations can be written as

σ_s(R,r) ≈ 2πRJ + 2πrε₋. (11)

Subtracting eqn (11) from eqn (10) we find the surface energy difference σ_b(R,r) − σ_s(R,r) ≈ −2rJ + πr(ε₊ − ε₋) which is independent of R. For symmetric impurity–solvent and impurity–solute interaction energy σ_b(R,r) − σ_s(R,r) ≈ −2rJ. This implies that the bonded cluster has less surface energy compared to separate clusters. Similar analysis has been carried out to calculate the free energy of a droplet doing cross-nucleation.³⁶

Fig. 10 A schematic diagram of cross nucleation to calculate the surface energy difference between the bonded and separate clusters.

Acknowledgements

We acknowledge the support from the EPSRC Programme Grant (Grant EP/R018820/1) which funded the Crystallization in the Real World consortium. In addition, we gratefully acknowledge the use of the computational facilities provided by the University of Warwick Scientific Computing Research Technology Platform.

Notes and references

1. R. Becker and W. Döring, Ann. Phys., 1935, 24, 752.
2. On the Theory of New Phase Formation. Cavitation, ed. R. A. Sunyaev, Princeton University Press, Princeton, 1992, pp. 120–137.
3. D. Kashchiev, Nucleation, Butterworth-Heinemann, Oxford, 2000.
4. S. Ryu and W. Cai, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011603.
5. D. Stauffer, A. Coniglio and D. W. Heermann, Phys. Rev. Lett., 1982, 49, 1299–1302.
6. V. A. Shneidman, K. A. Jackson and K. M. Beatty, J. Chem. Phys., 1999, 111, 6932–6941.
7. F. Schmitz, P. Virnau and K. Binder, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 053302.
8. N. Duff and B. Peters, J. Chem. Phys., 2009, 131, 184101.
9. Y. Lifanov, B. Vorselaars and D. Quigley, J. Chem. Phys., 2016, 145, 211912.
10. S. Teychené, I. Rodríguez-Ruiz and R. K. Ramamoorthy, Curr. Opin. Colloid Interface Sci., 2020, 46, 1–19.
11. D. A. Vorontsov, G. Sazaki, E. K. Titaeva, E. L. Kim, M. Bayer-Giraldi and Y. Furukawa, Cryst. Growth Des., 2018, 18, 2563–2571.
12. D. Gebauer, H. Cölfen, A. Verch and M. Antonietti, Adv. Mater., 2009, 21, 435–439.
13. S. Xu, Y. Liu, R. Zhang and Y. Wang, Cryst. Growth Des., 2024, 24, 103–114.
14. D. Han, Y. Wang, Y. Yang, T. Gong, Y. Chen and J. Gong, CrystEngComm, 2019, 21, 3576–3585.
15. J. Anwar, P. Boateng, R. Tamaki and S. Odedra, Angew. Chem., Int. Ed., 2009, 48, 1596–1600.
16. J. Anwar and D. Zahn, Angew. Chem., Int. Ed., 2011, 50, 1996–2013.
17. A. A. Bertolazzo, P. M. Naullage, B. Peters and V. Molinero, J. Phys. Chem. Lett., 2018, 9, 3224–3231.
18. G. G. Poon, T. Lemke, C. Peter, V. Molinero and B. Peters, J. Phys. Chem. Lett., 2017, 8(23), 5815–5820.
19. A. Hudait and V. Molinero, J. Am. Chem. Soc., 2014, 136, 8081–8093.
20. D. Mandal and D. Quigley, Soft Matter, 2021, 17, 8642–8650.
21. A. J. Page and R. P. Sear, Phys. Rev. Lett., 2006, 97, 065701.
22. L. O. Hedges and S. Whitelam, Soft Matter, 2012, 8, 8624–8635.
23. P. Grosfils and J. F. Lutsko, Crystals, 2021, 11, 4.
24. D. Mandal and D. Quigley, J. Chem. Phys., 2022, 157, 214501.
25. F. Ettori, T. J. Sluckin and P. Biscari, Phys. A, 2023, 611, 128426.
26. L. Yao and R. L. Jack, J. Chem. Phys., 2023, 159, 244110.
27. G. G. Poon, S. Seritan and B. Peters, Faraday Discuss., 2015, 179, 329–341.
28. G. Torrie and J. Valleau, J. Comput. Phys., 1977, 23, 187–199.
29. S. Auer and D. Frenkel, Annu. Rev. Phys. Chem., 2004, 55, 333–361.
30. R. J. Allen, C. Valeriani and P. R. ten Wolde, J. Phys.: Condens. Matter, 2009, 21, 463102.
31. F. A. Escobedo, E. E. Borrero and J. C. Araque, J. Phys.: Condens. Matter, 2009, 21, 333101.
32. R. J. Allen, P. B. Warren and P. R. ten Wolde, Phys. Rev. Lett., 2005, 94, 018104.
33. S. Ružička, D. Quigley and M. P. Allen, Phys. Chem. Chem. Phys., 2012, 14, 6044–6053.
34. K. E. Blow, G. A. Tribello, G. C. Sosso and D. Quigley, J. Chem. Phys., 2023, 158, 224102.
35. L. Onsager, Phys. Rev., 1944, 65, 117–149.
36. N. H. Fletcher, J. Chem. Phys., 1958, 29, 572–576.
37. L. Yu, J. Am. Chem. Soc., 2003, 125, 6380–6381.
38. S. F. S. P. Looijmans, D. Cavallo, L. Yu and G. W. M. Peters, Cryst. Growth Des., 2018, 18, 3921–3926.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00815d

---