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

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.


I. 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 4 .This model has been used to test assumptions and predictions of CNT in several previous studies [5][6][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 8 and more complex nucleation pathways 9 .
Impurities or additives are often used to control the nucleation process 10 , either accelerating 11 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 12 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 a) Electronic mail: dipanjan.mandal@warwick.ac.uk b) Electronic mail: d.quigley@warwick.ac.uk 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 17 .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 18 .Presence of ionic impurities in aqueous nano-droplet lead to demixing of low density amorphous ice and impurity rich aqueous glass when cooled 19 .
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 20 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 21 , pore-geometry 22 and surface roughness 23 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 24 provides insight into multi-component nucleation mechanism.Furthermore, examining the role of defects on magnetic droplet nucleation 25 , studying nucleation in the random field Ising model 26 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 27 , a surfactant-like additive was used to enhance the nucleation rate by decreasing the surface tension, similar to our previous work 20 .
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 II, 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 Section III and Section IV for static and dynamic impurities respectively.We characterize the different nucleation regimes that emerge from interaction with mobile impurities in Section IV.Section V 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 VI.Finally, we conclude in Section VII.

II. 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 as where ǫ 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 a solute reservoir and 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 (ǫ + = ǫ − ), antisymmetric (ǫ + = −ǫ − ) 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 ∆y = ± √ d 2 − ∆x 2 .Otherwise if ξ ≥ α we attempt to transmute a randomly selected solute or solvent.One Monte Carlo step consists of L 2 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][31][32][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 0 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 λ 0 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 0 which is system size dependent for a given choice of λ 0 , however the overall rate is not 34 .
The details of our implementation of both umbrella sampling and forward flux sampling algorithms could be found in Ref. 20 .

III. 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 solute-impurity interactions with an equal number of solvent-impurity 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 .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 B T 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 βF where B = − log ρ 1 − A s + A b , and ρ 1 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 Eq. 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 Eq. 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 impuri- 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 uppointing triangles.Similar behaviour is also observed for λc as a function of βǫ.Dashed horizontal lines represent the barrier height βF (λc) = 40.81and critical cluster size λc = 496 for the pure system (ρi = 0) at βJ = 0.67 and β∆µ = 0.067.
ties when βǫ d is large and positive.Here the free energy shows significant deviation from the expression shown in Eq. 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 βJ c = ln(1 + √ 2)/2 = 0.4406 . . . 35.
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 .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 Supporting Information (SI) for detailed free energy plots.
The barrier heights and critical sizes are consistent with our previous study when ǫ d = 0 20 .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 26 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.

IV. DYNAMIC IMPURITIES
In the case of dynamic impurities we include nonlocal impurity swap moves as described in section II.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.

A. 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).
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 SI 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 SI 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 27 .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 .

B. Symmetric interaction energy
Now we set the interaction energy ǫ + = ǫ − = ǫ to be symmetric.In this case the total solute-impurity and solvent-impurity interaction energy takes the constant value 4ǫρ i L 2 , 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 ρ i L 2 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 SI, we have plotted the variation of free energy for different negative values of symmetric interaction energy between impuritysolute 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 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.

C. Behaviour map for asymmetric interactions
If allowing the impurity-solute and impurity-solvent interaction energies to be neither symmetric nor antisymmetric, 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 nearestneighbour 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 20 .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.
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 V I), 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 (βǫ + + βǫ − , βǫ + − βǫ − ).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 .

V. 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 Eq. 3. The right hand side of Eq. 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.
The nucleation rates extracted from Fig. 7 are replotted 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.
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 subsection IV B, we do not see a significant variation in the nucleation rate for symmetric interactions where βǫ is negative and  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.
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 Eq. 4. In the case of dynamic impurities, we use A b and A s as fitting parameters to fit Eq.( 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 I.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 βǫ.
The Becker-Doring-Zeldovich expression of the nucleation rate can be written as where D c is the diffusion coefficient, Γ is the Zeldovich factor evaluated at the critical cluster size λ c .Using Eq. 4,  these quantities can be expressed as The diffusion coefficient D c = ∆λ(t) 2 /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 SI for D c estimation.In Table .II, we compare the nucleation rates obtained from the Becker-Doring-Zeldovich analysis (see Eq. 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 F F S , are the rates obtained from Eq. 5 and independent forward flux sampling simulations respectively.The maximum error in determining I BDZ is obtained by the expression where ∆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 Table .S1 and Table .S2 in the SI comparing I BDZ and I F F S at different regimes of the behaviour map.

VI. IMPURITY CLUSTERING AND CROSS-NUCLEATION
Impurities form multiple clusters when both impuritysolute 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.
Assuming a circular shape of both nuclei we may write the surface energy difference between the bondedcluster 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 36 .
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.

VII. 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 impuritysolvent 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.saturation criterion could be related with competition between the size of the largest void area without impurities and the critical cluster size.If the size of the largest void area is greater than the critical cluster size, we expect to see the saturation in the free energy barrier even for static impurities.Examining the detailed statistics of void site and distribution expected from a uniform distribution of impurities could in principle lead to an estimate of that threshold.For the impurity density ρ i = 0.004, we observe saturation both in free energy barrier [see Fig. S3(a)] and nucleation rate [see Fig. S4] with increasing βǫ d .However, at higher impurity density ρ i = 0.008, when the average void area excluding impurities decreases, we observe a monotonic increase in barrier height without saturation as βǫ d is increased [see Fig. S3(b)].We also note the differences in shape of saturated free energy between static and dynamic cases at ρ i = 0.004.Unlike dynamic impurities, the free energy curve becomes flatter in the case of static impurities and starts to deviate from the standard form of the free energy function assumed in classical nucleation theory as given in Eq.4 (see βǫ d = 1.6 curve in Fig. S7 and high positive βǫ d curves in Fig. S3).The confinement/constraint imposed by the immobile impurities could be responsible for this behaviour as it forces nuclei into shape with surface area to perimeter ratios that differ from the ideal case.
2. Free energy barrier with symmetric interaction energy (ǫ+ = ǫ− = ǫ): In the surfactant regime of the behaviour map we do not see much variation in the nucleation rate as shown in Fig. 8(a) of the paper, for dynamic impurities with βJ = 0.83, β∆µ = 0.083, ρ i = 0.012 and α = 0.05.Similar behaviour is reflected in free energy plots for different values of symmetric interaction energies that belong to the surfactant regime as shown in Fig. S5.

Decrement in free energy barrier height due to mobile impurities:
As observed in Section 4, the free energy barrier height to nucleation decreases when impurities are dynamic compared to the static case for same set of interaction energies as shown in Fig. S6 (see Fig. 2(a) of the paper).In this case the interaction energies are anti-symmetric (βǫ + = −βǫ − = 0.4) and other parameter values are βJ = 0.67, β∆µ = 0.067 and ρ i = 0.008.Dynamic impurities enhance the nucleation rate by decreasing the barrier height.In this example the microscopic interaction between impurity-solute and impurity-solvent are respectively weakly-repulsive and weakly-attractive and lies at the boundary of the surfactant regime of the behaviour map.We also observed similar decrement in barrier height in our earlier work [D.Mandal and D. Quigley, Soft Matter, 2021, 17, 8642-8650] for impurities with neutral interactions.However, in that case, because of the neutral interaction impurities prefer to stay at the boundaries of the cluster, it reduces the interfacial free energy.S3.Nucleation free energy with varying dimensionless interaction energy difference βǫ d with fixed βJ = 0.67, β∆µ = 0.067 for system size L = 100 at static impurity density (a) ρi = 0.004 and (b) ρi = 0.008.We see no further increase in free energy barrier height with increasing βǫ d when ρi = 0.004 or lower, i.e., when the impurities are sparsely distributed so that a critical cluster can fit in the void space without interacting with impurities.This behaviour in barrier height is not observed for ρi = 0.008 when impurity density is higher.Free energy barrier for the system without impurities (ρi=0) is plotted by dotted line for comparison.

Fitting free energy barrier:
Fitting of free energy to the expression given in Eq. 4 in the case of static impurities with density ρ i = 0.004, βJ = 0.67 and β∆µ = 0.067 for different βǫ d is shown in Fig. S7, where we allow the surface A s and bulk A b terms to vary from the ρ i = 0 case.We see that A s increases and A b decreases monotonically with increasing βǫ d from negative to positive values.The fitting becomes more accurate with decreasing βǫ d .
Similar fitting of calculated free energy with Eq. 4 for dynamic impurities with βJ = 0.67, β∆µ = 0.067, ρ i = 0.004 and α = 0.05 is shown in Fig. S8 for different βǫ.We see monotonic decrease and monotonic increase of the bulk and surface terms respectively with increasing βǫ until they converge to non-zero finite values after entering into the inert-spectator regime of the behaviour map.We note that, for pure Ising model at low temperatures ∆g ≈ ∆µ, and from Fig. S8 we see that this relation holds for neutral impurities, as A b ≈ β∆µ when βǫ = 0, but not for non-zero interaction energies.

Nucleation rate and diffusion coefficient Dc:
The Becker-Doring-Zeldovitch nucleation rate I BDZ is calculated using Eq. 5 of the paper for high ρ i = 0.02 and intermediate ρ i = 0.012 impurity densities.For that, the diffusion coefficient D c is calculated after performing independent simulations starting from the critical cluster size at time t = 0 and calculating the slope of the mean squared deviation of cluster size with time t.Estimated values of D c corresponding to plots displayed in Fig.    are written in column 4 of Table S1.The values of different parameters required for calculating I BDZ is given in Table S1 and Table S2 for ρ i = 0.02 (anti-symmetric interaction energy) and ρ i = 0.012 (symmetric interaction energy)   ) for βJ = 0.83, β∆µ = 0.083, ρi = 0.012 and α = 0.05 with symmetric interaction energy βǫ+ = βǫ− = ǫ.This range of βǫ belongs to the regime in which impurities act like heterogeneous nucleating sites.
respectively.In the final column the nucleation rate obtained from independent forward flux sampling simulations I F F S is written.The results for I BDZ and I F F S match quite well for the range of interaction energies considered in both tables.

FIG. 1 .
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.

FIG. 4 .
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.

11 FIG. 5 .
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.

FIG. 8 .
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 Is(ρ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.

10
FIG. S4.Rate of obtaining a cluster of size λ starting from a metastable solution phase for different interaction energy difference βǫ d at fixed βJ = 0.67, β∆µ = 0.067 and ρi = 0.004 with static impurities.The constant value of I(λ) for large λ is the nucleation rate.
FIG. S9.Linear fitting of diffusivity ∆λ 2 = 2Dct obtained from simulation for dynamic impurities with βJ = 0.83, β∆µ = 0.083, ρi = 0.02 and α = 0.05.Initial size of the cluster is set to the critical cluster size calculated from the position of the maxima in the respective F (λ) vs. λ plots.

TABLE S2 .
Comparison of nucleation rates obtained from Becker-Doring-Zeldovich analysis (IBDZ) and forward flux sampling method (IF F S