Nucleation rate in the two dimensional Ising model in the presence of random impurities

Abstract

Nucleation phenomena are ubiquitous in nature and the presence of impurities in every real and experimental system is unavoidable. Yet numerical studies of nucleation are nearly always conducted for entirely pure systems. We have studied the behaviour of the droplet free energy in two dimensional Ising model in the presence of randomly positioned static and dynamic impurities. We have shown that both the free energy barrier height and critical nucleus size monotonically decreases with increasing the impurity density for the static case. We have compared the nucleation rates obtained from the Classical Nucleation Theory and the Forward Flux Sampling method for different densities of the static impurities. The results show good agreement. In the case of dynamic impurities, we observe preferential occupancy of the impurities at the boundary positions of the nucleus when the temperature is low. This further boosts enhancement of the nucleation rate due to lowering of the effective interfacial free energy.

---

1 Introduction

Nucleation phenomena are of great importance both in physics and chemistry. The study of nucleation has a long history, rooted in the widely used theoretical model of classical nucleation theory (CNT)¹ and the Becker–Döring (BD) expression for the nucleation rate.² CNT is able to explain the qualitative (and sometimes quantitative) behaviour, of various nucleation processes occurring in natural and experimental systems.^3,4 The nucleation rate at which a stable phase nucleates from a metastable parent depends on various parameters, e.g., temperature,⁵ pressure,⁶ external electric field,⁷ viscosity,⁸ application of shear,⁹etc. It is also heavily influenced by surfaces and the presence of impurities^10,11 in the system.

Modelling nucleation from solution is an active area of research with relevance to pharmaceutical manufacture, biomineral formation and other highly complex precipitation processes. Atomistic simulations¹² are now able to make quantitative predictions of the nucleation rate. Unfortunately, even for nucleation of simple salt crystals from solution^13,14 agreement with experiment remains elusive and sensitive to details of the atomistic force field and the order parameter used to identify solid regions in the simulation.¹⁵ In addition, the solution is nearly always simulated as pure with no impurities or spectator/counter ions. A quantitative understanding of how these additional species impact on the nucleation rate is lacking.

The most basic physics of nucleation from solution can be captured using a simple lattice-based model of solute precipitation.¹⁶ Minimal models of this kind in two and three dimensions have been used to probe assumptions of CNT and to test advanced methods for quantifying rare events against well-established results.^17–21 Extensions of the basic Ising-like model have been used to study nucleation inside a rectangular pore,^22,23 on substrates,²⁴ two-step nucleation mechanisms,²⁵ and nucleation in the presence of shear⁹ and on single immobile microscopic impurities.²⁶ Studies have also explored the choice of microscopic kinetics on the nucleation rate.²⁷ In this work we study such a model in the presence of neutral impurities as a further step toward introducing some of the complexity found in “real” solutions. We explore how this changes nucleation rate, critical nucleus size and free energy barrier height, as a function of impurity concentration.

We consider two scenarios, static and dynamic impurities. In the static impurity case the impurities are sites distributed randomly in space and fixed throughout the simulation. We work with quantities averaged over many samples of this static disorder. In the dynamic case impurities can migrate via Kawasaki dynamics.²⁸ A particular advantage of working with such simplified models is that the timescale of impurity migration relative to cluster growth can be adjusted to probe different kinetic regimes.

We have used umbrella sampling techniques to calculate the free energy barrier for different densities of the impurity sites. We also calculate the nucleation rate using the Forward flux sampling (FFS) technique and compare the results to those computed from the free energy barrier via the Becker–Döring expression. We have shown that the impurity particles preferentially occupy the boundary positions of the nucleus for the dynamic case when the temperature is low. This surface accumulation process enhances the nucleation rate by a multiple of 10⁴ compared to the static case.

The rest of this paper is organised as follows. In Section 2, we describe the model and simulation techniques used to study the model. Variation of the free energy barrier height and the critical nucleus size are shown in Section 3. Section 4 contains the comparison of theoretical and computational nucleation rates for different impurity concentration. The effect of dynamic impurities on barrier height as well as nucleation rates are described in Section 5. Finally we conclude in Section 6.

2 Model & simulation techniques

We consider the two dimensional Ising model on a L × L square lattice with spin S_i at position i. The Hamiltonian of the model may be written aswhere J is the strength of the coupling between two nearest neighbours, h is the externally applied field and 〈i,j〉 represents summation over all nearest-neighbour pairs in the usual manner. The spin can take two values, S_i = ±1. We take the strength of the coupling to be positive (J > 0), i.e., the ferromagnetic Ising model. In our simulations we have taken the copling strength J = 1. This is analogous to a lattice-gas model in which each lattice site can be occupied (S_i = +1) or unoccupied (S_i = −1), with the external field playing the role of a particle reservoir at a constant chemical potential.

We study the discontinuous magnetisation reversal for a system with initially negative magnetisation in the presence of a positive magnetic field. This necessarily occurs only for temperatures T < T_c (T_c is the critical temperature) and proceeds via nucleation and growth of domains with positive magnetisation. In the lattice-gas (solute precipitation) analogy this corresponds to an initially mostly-empty lattice (solvent-rich) to a full lattice (solute-rich) at a chemical potential where equilibrium with a particle reservoir corresponds to a mostly-occupied lattice.

We modify the model by introducing a third spin state S_i = 0 at random positions of the lattice. This represents non-magnetic impurities, or “neutral” impurity ions in solution which do not favour interaction with either solute or solvent neighbours. The interaction energy of impurity particles with both solute and solvent is chosen to be half way between solute–solvent and solute–solute (or solvent–solvent) interaction energy, i.e. zero. Other possibilities could be studied by creating impurities which interact more favourably with solvent than solute (for example). In the present work we restrict ourselves to the neutral case. These non-zero spins are evolved in time via the usual spin–flip dynamics in which a randomly selected spin is flipped, and the move accepted or rejected according to the Metropolis criterion. In what follows we mainly use the magnetic terminology for brevity.

The overall density of spin-0 impurities is kept fixed throughout the simulation. Further, we consider two different circumstances, static and dynamic impurities. In the static case, the impurity particles are immobile. In the dynamic case impurities can diffuse through the lattice using spin-exchange dynamics often known as Kawasaki dynamics.²⁸ The ratio of impurity-diffusion move attempts to spin flip move attempts controls the mobility of impurities relative to the nucleation timescale.

We calculate the nucleation rate using two independent methods, and as a function of the impurity density. The first method calculates the free energy barrier to nucleation and applies the BD expression for the nucleation rate and the second method is the Forward Flux Sampling (FFS) method which calculates the rate independently of classical nucleation theory. A typical configuration of the system during the transition is shown in Fig. 1 where red represents occupied (S_i = 1) sites, blue represents impurity (S_i = 0) and white represents empty sites (S_i = −1).

Fig. 1 Typical configuration of the nucleating system in the presence of static impurities for T = 1.5, h = 0.05, ρ_i = 0.028 and L = 100. Particles and impurities are respectively denoted by red and blue colors respectively, and white represents empty site.

In the next sections, we describe in more detail the simulation techniques used to calculate the barrier heights and the nucleation rates.

2.1 Umbrella sampling

In situations where the free energy barrier to nucleation is large, it is impossible to adequately sample the probability distribution for large nuclei in an unbiased simulation. Umbrella sampling (US) is widely used biased simulation technique to calculate free energy as a function of a collective variable or reaction coordinate. The US method was introduced by Torrie and Valleau in 1976²⁹ and tested successfully for Lennard-Jones systems. The method has been used in different contexts, e.g., protein folding dynamics in biological systems, measuring reaction rates in chemical systems, etc.³⁰ The method is also often used (as here) in the context of nucleation^31,32 to obtain free energy as a function of cluster size and hence parameterise a BD calculation of the nucleation rate. We define a cluster as a contiguously connected region of +1 spins. This definition has proven to yield accurate BD estimates of the nucleation rate in the regime of interest here,¹⁸ but may not be optimal closer to T_c.³³ The configuration space is divided into overlapping windows of equal size which span the range of cluster sizes from one to the nucleus size greater than the critical nucleus size. The probability distribution of the cluster sizes lying inside the window range is calculated for each window. In our US simulations all windows span a cluster size of range 20. A range of 10 overlaps between two neighbouring windows. We have taken wide overlap between two windows to reduce the error while combining different parts of the free energies. We use an infinite square well umbrella potential of the same width as the window. For each window we run the simulation for 10⁹ steps where each step is an attempt to flip a randomly chosen spin. We skip the spin update if the the random site contain an impurity spin. After each step we measure cluster sizes. If the largest cluster size escapes the window range we reject the spin flip and restore the previous state of the spin. Simulations of this kind are particularly suited to parallelisation over windows, for which we employ OpenMP³⁴ threading to decrease simulation time significantly. In the case of static impurities, we average the obtained probability distribution for each window with 28 different impurity configurations. An uncertainty in the probability distribution is obtained by averaging over 28 independent calculations of probability distribution in each window. The maximum error is obtained from the standard error of the probability at each nucleus size, where σ is the standard deviation of [scr N, script letter N] independent results.

The relative free energy for a droplet of size λ at the n-th window may be written as

f^US_n(λ) = −k_BT ln[P_n(λ)], (2)

where P_n(λ) is the probability of sampling a cluster of size λ and k_B is the Boltzmann constant. We take k_B = 1 in our calculations. We combine the free energy obtained for the n-th window with the free energies of the previous windows by a constant shift S_n which may be written aswhere m is the total number of overlapped points and i runs over all overlapped points between (n − 1)-th and n-th windows. After the constant shift the free energy of the n-th window may be written as

F^US_n(λ) = f^US_n(λ) + S_n. (4)

Once all windows have been combined in this fashion we obtain F^US(λ), the free energy over the entire sampled range of cluster sizes.

2.2 Forward flux sampling

FFS^35–38 is a numerical technique widely used to calculate the rate of occurrence of a rare event. In this paper we calculate the nucleation rate from the initial metastable phase (negative spin) to the stable phase (positive spin) using FFS. The method proceeds via the following steps. First, we divide the configuration space connecting the metastable and the stable phase with intermediate equally spaced interfaces λ_i. Each interface is an isosurface of the maximum cluster size, such that progression from the first interface λ₀ to the final interface λ_N captures a nucleation trajectory. We assume that λ_i+1 is always greater then λ_i.

The nucleation rate is then estimated as

where I₀ is the initial positive flux through interface λ₀, i.e., the number of crossings of λ₀ per unit area in unit time. We define a unit of time as L × L attempts to flip a single spin. The units of I₀ and hence I^FFS are therefore crossings per Monte Carlo step per single site, or MCSS⁻¹. P(λ_i+1|λ_i) is the probability that a sequence of Monte Carlo moves beginning with a configuration with maximum cluster size λ_i will cross to λ_i+1 before reaching λ_A = 8, i.e., reaches the next interface before returning to the metastable parent phase (λ ≤ λ_A).

In our simulations, we have used 250 interfaces with a constant gap of 10 between adjacent interfaces. The position of the first interface varies for different temperatures. For T = 1.5 and h = 0.05, the first interface position is λ₀ = 16. The position of first interface is chosen such that crossings are infrequent, i.e. the probability of being at λ = λ₀ state is quite low (λ₀ is above the 95th percentile of the cluster probability in the parent phase). We run trials starting from the metastable (negative spin) phase to calculate the initial positive flux, storing configurations which cross λ₀ in the direction of increasing λ. We then initialise Monte Carlo trajectories from randomly chosen configurations at this interface and count how many reach λ₁ before returning to λ_A to obtain the crossing probability of the first interface. This is repeated for all subsequent interface pairs. In each case we sample the interface until 25 200 successful events are captured. The initial flux is calculated based on how much time is needed to produce 25 200 positive crossings of λ₀. In a typical FFS simulation we sample interface up to λ ∼ 2500. Plots of the interface probabilities at λ = λ_i and the rate of reaching state λ starting from the metastable phase for different impurity densities ρ_i are shown in Fig. 2(a) and (b) respectively for fixed T = 1.5 and h = 0.05. The interface probability and rate saturate for large λ indicating post-critical behaviour. The saturation becomes slower with increasing ρ_i.

Fig. 2 Plots of (a) interface probability P(λ_i+1|λ_i) and (b) nucleation rate I_λ to a cluster of size λ as a function of λ for different impurity densities ρ_i with fix temperature T = 1.5 and field h = 0.05.

3 Barrier height and the critical nucleus size

We first report on variation of the nucleation behaviour in the presence of randomly distributed static impurity spins.

The free energy of a cluster of +1 spins of size λ may be written as

where ρ₁ is the monomer density, i.e., the fraction of the sites occupied by isolated positive spins in the metastable parent phase and P(λ) is the probability that an observed cluster will be of size λ. For a cluster size upper limit of λ = λ_m, this probability obeys the normalisation conditionThe variation of the free energy F^US(λ) for different densities of the impurities ρ_i is shown in Fig. 3. The height of the free energy barrier decreases monotonically with increasing ρ_i. The nucleus size at which the barrier height takes maximum value is known as the critical nucleus size λ_c which monotonically decreases with increasing ρ_i. This suggests that the presence of spin-0 impurities enhances nucleation of S_i = +1 domains.

Fig. 3 Plot of the free energy F^US(λ,ρ_i) obtained from the US simulations as a function of the nucleus size λ for different static impurity densities ρ_i at T = 1.5 and h = 0.05. The black line is the fit to the modified CNT expression of the free energy F(λ,ρ_i) (see eqn (9)). The error-bars are less than point size.

To gain insight into this enhancement, we fit our free energy barriers to nucleation theory and extract trends in the resulting physical parameters. In a previous study Ryu and Cai¹⁸ have numerically calculated free energy barriers in the absence of impurities and confirmed these accurately obey a cluster free energy given by

where A₁ is proportional to the interfacial tension per unit length, A₂ is the coefficient of the logarithmic correction term introduced by Langer³⁹ to incorporate the microscopic fluctuation of the cluster shape. In ref. 18, the term A₃ is obtained from matching F(λ) with the analytical expression of droplet free energy⁴⁰ for small λ. However, in the presence of impurities, we have used modified expression of free energy which can be written as

A₃(ρ_i) = −k_BT ln ρ₁ − A₁(ρ_i) + 2h. (10)

The expression of A₃(ρ_i) is obtained from the equality F(1,ρ_i) = −k_BT ln ρ₁, where ρ₁ is the density of isolated +1 spins in the presence of impurities. The term A₂ can be calculated theoretically for homogeneous nucleation and written as .⁴¹ We fit the free energy plot of the system without impurities obtained from the umbrella sampling techniques to the function F(λ,ρ_i) written in eqn (9) to estimate A₁(ρ_i) and A₃(ρ_i). The fitted values of the parameters A₁(ρ_i) and A₃(ρ_i) are respectively 4.279 and 3.776 for T = 1.5 and h = 0.05, which are in very close agreement with A₁ ≈ 4.3 and A₃ ≈ 3.7, values obtained in ref. 18 for ρ_i = 0.

In the presence of impurities we fix for all ρ_i. Allowing this parameter to vary has negligible impact on the quality of fit suggesting that the low impurity densities we study here do not significantly alter fluctuations in cluster shape.

The size-independent term A₃(ρ_i) is obtained from the density of isolated +1 spins during unbiased simulations of the system in the metastable phase. This varies with impurity concentration. The only free parameter when fitting the free energy plots in Fig. 3 for different ρ_i is the interfacial tension per unit length A₁(ρ_i). As shown in Fig. 4(a), this decreases linearly with increasing impurity density for fixed T and hence can be written as

A₁(ρ_i) = m₁ρ_i + m₂. (11)

It is interesting to explore how the parameters m₁ and m₂ vary with temperature and field. We observe that A₁(ρ_i) increases with decreasing T for fixed ρ_i. The variation of A₃(ρ_i) as a function of ρ_i for T = 1.4 and 1.5 is shown in Fig. 4(b). A₃(ρ_i) increases with increasing ρ_i for both temperatures.

Fig. 4 Variation of (a) the interfacial tension A₁(ρ_i), and (b) the size-independent term A₃(ρ_i) as a function of the impurity densities ρ_i for T = 1.5 and h = 0.05. The maximum impurity density is ρ_i = 0.028 i.e. 2.8% of the lattice is occupied by the impurities.

Spin-0 impurities which bridge between sites occupied by opposite spins lower the overall energy of the system by a greater amount than impurities surrounded entirely by spins of the same sign. Hence it is the presence of spin-0 impurities at the cluster boundary (rather than inside the cluster) which has the greatest impact on the nucleation barrier. With static impurities, the average impurity density per unit length of cluster circumference is independent of both temperature and cluster size, and hence we can capture their impact on the interfacial free energy via this simple functional form.

Now we compare the critical nucleus size and the barrier height obtained from our modified free energy expression and from the simulation. From eqn (8) it is easy to show that the critical nucleus size may be written as

The variation of λ_c and F(λ_c,ρ_i) for different ρ_i obtained from theory and simulations are shown in Fig. 5(a) and (b) respectively for two different T. The results agree satisfactorily. As would be expected from the reduced interfacial free energy, increased density of impurities results in a smaller critical nucleus size and barrier height.

Fig. 5 Variation of (a) the critical nucleus size λ_c and (b) the maximum free energy barrier F(λ_c,ρ_i) as a function of ρ_i for T = 1.5 and h = 0.05.

4 Nucleation rates

We have calculated nucleation rates using both BD theory, and the FFS technique (which is independent of CNT assumptions). The expression of the nucleation rate using FFS is shown in eqn (5). We compare these nucleation rates obtained from FFS with nucleation rate estimates from Becker–Döring theory using our modified expression for F(λ_c,ρ_i) (see eqn (9) and (12)) which incorporates the effects of impurities. The nucleation rate may be written aswhere D_c is a diffusion coefficient which captures the rate of addition to the critical nucleus and Γ is the Zeldovich factor defined as

The quantity D_c is extracted from the mean squared variation in the nucleus size as plotted in Fig. 6(a) with increasing time for different impurities densities. We perform separate Monte Carlo simulations to calculate D_c for different ρ_i. The initial configuration at t = 0 is taken to be the nucleus with critical size. The mean squared variation is defined as

Δλ²(t) = [λ(t) − λ_c]². (15)

Fig. 6 Variation of 〈Δλ²(t)〉 as a function of time t with (a) static impurities for different densities and (b) dynamic impurities for different α.

We average Δλ²(t) over an ensemble of 10⁵ such trajectories for each line plotted in Fig. 6(a). The diffusion coefficient D_c is defined as

where 〈〉 represents the ensemble average. The slope gradually decreases with increasing ρ_i indicating that the presence of impurities hinders addition of new positive spins to the cluster. Since this term appears in pre-factor of the exponential it has little impact on the nucleation rate, such that the reduction in F(λ_c,ρ_i) dominates leading to an overall increase in nucleation rate with increasing impurity density.

The nucleation rates obtained from the FFS [see eqn (5)] and BD theory [see eqn (13)] are compared in Fig. 7 which shows good agreement. From this we conclude that the impact of static impurities on nucleation rate is well described via classical nucleation theory with only minor extensions required to capture trends with increasing impurity density over the range explored 0 ≤ ρ_i ≤ 0.028. For higher impurity densities we observe multiple overlapping pre-critical nuclei and therefore CNT is not appropriate.

Fig. 7 Comparison of the nucleation rates obtained from the FFS method and the Becker–Döoring theory for T = 1.5 and T = 1.4 with h = 0.05.

We have plotted the nucleation rate I(MCSS⁻¹) as a function of ρ_i for different T with fixed h = 0.05 in Fig. 8(a). The nucleation rate grows exponentially with increasing impurity density and the exponential growth becomes faster with decreasing temperature. The slope of the best fitted line increases with decreasing T for fixed h as shown in Fig. 8(a). This indicates that the barrier height decreases more rapidly with increasing ρ_i when the temperature is lowered. The variation of nucleation rates for different h with fixed T = 1.5 is shown in Fig. 8(b). The rates grow exponentially as a function of ρ_i for fixed T and become more steeper with decreasing h. We can approximately write the nucleation rate as I ≈ C exp(ν_Tρ_i) and I ≈ C exp(ν_hρ_i) for fixed h and T respectively. Both ν_T and ν_h increase monotonically with decreasing T and h.

Fig. 8 Variation of the nucleation rate as a function of ρ_i for different values of the temperatures with fixed external field h = 0.05.

5 Dynamic impurities

With static spin-0 impurities, the reduction of the interfacial free energy is controlled by the overall impurity density. The probability of a impurity being located at the growing +1 spin boundary is constant, and independent of temperature. With dynamic impurities, the possibility of two limiting regimes arises.

(1) Low mobility impurities: the boundary of the growing +1 spin cluster will expand rapidly compared to the randomly distributed and slow-moving impurities, resulting in a reduction of the interfacial free energy similar to the static case.

(2) High mobility impurities: sufficiently fast-moving impurities can reach a local equilibrium distribution much faster than the +1 spin nucleus can grow or shrink. This implies that a description of the nucleation process in terms of a free energy at each nucleus size remains appropriate.

Here we apply both the BD and FFS approaches in the high mobility regime and compare to the static regime studied above.

We consider dynamic impurities of fixed density ρ_i. Impurity sites may move to the one of the four neighbouring sites via Kawasaki dynamics. The spin 0 impurity (blue) may exchange with a spin +1 (red) or spin −1 (white) neighbouring sites. We do not consider spin exchange dynamics involving only red and white sites. We define the quantity α such that the mobility of impurity particles increases with increasing α. The mobility parameter α takes values between 0 to 1. In our Monte Carlo simulations we attempt a spin exchange move for a random impurity particle with probability α. Attempts to flip the spin of randomly chosen size with spin ±1 are attempted with probability 1 − α.

Rates are reported to be consistent with the previous section, i.e., based on a unit of time during which, on average, every site on the lattice is subjected to one trial flip. Hence for large α, many Kawasaki exchange moves involving impurities are attempted per time unit.

We have plotted the free energy barrier to nucleation both for static (α = 0) and high mobility (α = 0.95) impurities with density ρ_i = 0.020 for different values of α. The free energy plots for T = 1.0 and T = 1.5 are shown in Fig. 9(a) and (b) respectively. We observe that the barrier height decreases significantly for the nonzero value of α compared to the α = 0 case when the temperature is low (T = 1.0). However, the critical nucleus size is approximately the same. By contrast, we see that the barrier height and critical nucleus size of the free energy do not change as dramatically between α = 0 and α = 0.95 at the higher temperature T = 1.5.

Fig. 9 Free energy barrier to nucleation with both static (α = 0) and high mobility (α = 0.95) impurities at (a) T = 1.0 and (b) T = 1.5 with h = 0.05 and ρ_i = 0.020. The critical nucleus size at T = 1.0 for α = 0.0 and α = 0.95 are respectively 756 and 665.

Nucleation rates I as a function of ρ_i obtained using FFS for static and high mobility impurities at T = 1.5 and h = 0.05 are shown in Fig. 10. We see that the rates obtained using FFS for different α are similar (i.e. within one order of magnitude) for every ρ_i at this high temperature, but there is an indication of high mobility impurities increasing the nucleation rate. This is consistent with the observation of marginally lower barrier heights for high α at T = 1.5 in Fig. 9(b).

Fig. 10 Comparison of nucleation rates obtained via FFS for static (α = 0) and high mobility (α = 0.95) impurities.

Based on the free energies in Fig. 9(a), we expect a much more dramatic impact of high mobility impurities at the lower temperature of T = 1.0. The FFS method becomes extremely inefficient at this temperature due to a very low density of spin +1 sites in the metastable phase. Computing a sufficiently accurate flux through λ₀ is intractable in this case.

We have, however, calculated the nucleation rate at T = 1.0 using the Becker–Döring expression. The rates are 4.005 × 10⁻³⁸ MCSS⁻¹ and 2.180 × 10⁻³⁴ MCSS⁻¹ for α = 0 and α = 0.95 respectively. The mobility of impurities hence enhances the nucleation rate by several orders of magnitude. The mean squared deviation 〈Δλ²(t)〉 for dynamic impurities with different α at T = 1.0 is shown in Fig. 6(b). We see that the slope of the plot reduces for α ≠ 0 implying the reduction of diffusion coefficient D_c for dynamic impurities.

The BD rates at T = 1.5 are 1.128 × 10⁻¹⁵ MCSS⁻¹ and 1.930 × 10⁻¹⁵ MCSS⁻¹ for α = 0.95 and α = 0.0 respectively. The corresponding critical cluster sizes are λ_c = 368 and λ_c = 363 for dynamic and static case respectively.

The role of high mobility impurities at high versus low temperature for three different temperatures with fixed h = 0.05 is illustrated in Fig. 11(a–f). Fig. 11(a, b), (c, d) and (e, f) are the snapshots of typical configuration after equilibrium at T = 1.5, T = 1.0 and T = 0.9 respectively for both static and dynamic impurities. For the purposes of this illustration we set the mobility parameter α = 0.95 prioritising the spin exchange dynamics for the impurity particles. Initially the impurities are placed randomly on the lattice. We create a circular nucleus of red (spin +1) particles at the centre of the lattice of size λ = 600 and evolve the system such that the cluster size remain confined within a range [λ − 10, λ + 10]. We see that the impurity particles preferentially occupy the boundary positions of the nucleus at low temperatures. This happens because the particular impurity particles modelled here act as a surfactant, reducing the surface free energy. However, at high temperatures the phenomenon of preferential occupancy of the blue particles at the boundary positions is entropically unfavourable leading to a less dramatic enhancement of the nucleation rate.

Fig. 11 Snapshot of the system after equilibration with static (left column) and dynamic impurities (right column) at (a and b) T = 1.5, (c and d) T = 1.0 and (e and f) T = 0.9 with fixed h = 0.05 and ρ_i = 0.020. The value of α is 0.95 for the dynamic impurities.

To quantify this observation, we plot the equilibrium fraction ϕ of the impurity particles present at the boundary of a static nucleus as a function of T for different sizes of the nucleus as shown in Fig. 12. The quantity ϕ is defined as

where N_r is number of red particles connected to the nucleus through the nearest neighbour sites and N_b is number of blue particles present at the boundary of the nucleus and also connected to the nucleus through the nearest neighbour sites. We set the mobility parameter α = 1 such that only spin exchange dynamics are performed for the impurity particles. The boundary particles are identified by measuring the distance r of the blue particle from the centre of mass of the nucleus. If r ≥ 0.7R, we count the blue particle as boundary blue particle, where R is mean radius of the nucleus. The quantity ϕ increases with decreasing T indicating the adsorption phenomenon of impurity particles at the boundary of the nucleus. This adsorption phenomenon is absent at high temperatures.

Fig. 12 Fraction of the blue particles present at the boundary of the nucleus with increasing the temperature for different size of the initial nucleus λ_i.

6 Conclusion

We have studied nucleation rates in the two dimensional Ising model in the presence of randomly placed impurities on a square lattice using two independent methods which are BD theory and FFS. The rates obtained using both methods agree satisfactorily for different impurity densities up to ρ_i = 0.028, i.e., when the only 2.8% of the sites are occupied by the impurities. The nucleation rate increases with increasing the impurity density. We have also studied the effect on the nucleation rates when the impurity particles are dynamic. The dynamics of the impurity particles accelerate the nucleation when the temperature is low. However, the impurity dynamics does not play as significant a role when the temperature is high.

We have not explored the regime of intermediate impurity mobility. Here one would not expect the distribution of impurities to equilibrium on a timescale more rapid than that on which spin +1 cluster grow or shrink. This would invalidate a free-energy based description, but could be amenable to a description based on BD theory if growth/shrinkage rates could be evaluated directly. Rates in this regime can possibly be probed via FFS, however in the high temperature regime accessible to this method we expect little impact of mobility on nucleation rate due to the less-dramatic difference between the two limiting cases of zero and high mobility.

We stress that our study has only considered the case where the impurities interact with both spin up and spin down sites in the same way. Different behaviour would be expected if (for example) the interactions favoured location of impurities within the bulk of the metastable phase. Future work will attempt to map behaviour in different interaction regimes such that approximate connections to solute precipitation can be made.

The study of the nucleation in the presence of random impurities on triangular lattice and three dimensional cubic lattice is an open area for future research. It would be interesting to investigate how the different dynamics of the impurity particles in other lattice geometries impact the nucleation rate.

Author contributions

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

Data availability

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

Conflicts of interest

The authors have no conficts of interest to declate.

Acknowledgements

We thank I. J. Ford for helpful discussions. We acknowledge the support from the EPSRC Programme Grant (Grant EP/R018820/1) which funds 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.

References

1. D. Kashchiev, Nucleation, Butterworth-Heinemann, Oxford, 2000.
2. R. Becker and W. Döring, Ann. Phys., 1935, 24, 752.
3. E. Mendez-Villuendas and R. K. Bowles, Phys. Rev. Lett., 2007, 98, 185503.
4. P. H. Poole, F. Sciortino, U. Essmann and H. E. Stanley, Nature, 1992, 360, 324–328.
5. A. Pitto-Barry and N. P. E. Barry, Angew. Chem., Int. Ed., 2019, 58, 18482–18486.
6. C. P. Lee and T. G. Wang, J. Appl. Phys., 1992, 71, 5721–5723.
7. D. Kashchiev, J. Cryst. Growth, 1972, 13–14, 128–130.
8. X. Frank, N. Dietrich, J. Wu, R. Barraud and H. Z. Li, Chem. Eng. Sci., 2007, 62, 7090–7097.
9. R. J. Allen, C. Valeriani, S. Tanase-Nicola, P. R. ten Wolde and D. Frenkel, J. Chem. Phys., 2008, 129, 134704.
10. R. M. Ginde and A. S. Myerson, J. Cryst. Growth, 1993, 126, 216–222.
11. L. Keshavarz, R. R. E. Steendam, M. A. R. Blijlevens, M. Pishnamazi and P. J. Frawley, Cryst. Growth Des., 2019, 19, 4193–4201.
12. V. Agarwal and B. Peters, Advances in Chemical Physics, John Wiley & Sons, Ltd, 2014, vol. 155, pp. 97–160.
13. H. Jiang, A. Haji-Akbari, P. G. Debenedetti and A. Z. Panagiotopoulos, J. Chem. Phys., 2018, 148, 044505.
14. N. E. R. Zimmermann, B. Vorselaars, D. Quigley and B. Peters, J. Am. Chem. Soc., 2015, 137, 13352–13361.
15. N. E. R. Zimmermann, B. Vorselaars, J. R. Espinosa, D. Quigley, W. R. Smith, E. Sanz, C. Vega and B. Peters, J. Chem. Phys., 2018, 148, 222838.
16. K. Binder and P. Virnau, J. Chem. Phys., 2016, 145, 211701.
17. H. Vehkamäki and I. J. Ford, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 6483–6488.
18. S. Ryu and W. Cai, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011603.
19. P. A. Rikvold, H. Tomita, S. Miyashita and S. W. Sides, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1994, 49, 5080–5090.
20. Y. Lifanov, B. Vorselaars and D. Quigley, J. Chem. Phys., 2016, 145, 211912.
21. D. Stauffer, A. Coniglio and D. W. Heermann, Phys. Rev. Lett., 1982, 49, 1299–1302.
22. A. J. Page and R. P. Sear, Phys. Rev. Lett., 2006, 97, 065701.
23. L. O. Hedges and S. Whitelam, Soft Matter, 2012, 8, 8624–8635.
24. M. L. Trobo, E. V. Albano and K. Binder, J. Chem. Phys., 2018, 148, 114701.
25. N. Duff and B. Peters, J. Chem. Phys., 2009, 131, 184101.
26. R. P. Sear, J. Phys. Chem. B, 2006, 110, 4985–4989.
27. J. Kuipers and G. T. Barkema, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011128.
28. K. Kawasaki, Phys. Rev., 1966, 145, 224–230.
29. G. Torrie and J. Valleau, J. Comput. Phys., 1977, 23, 187–199.
30. C. Matthews, J. Weare, A. Kravtsov and E. Jennings, Mon. Not. R. Astron. Soc., 2018, 480, 4069–4079.
31. S. Auer and D. Frenkel, Annu. Rev. Phys. Chem., 2004, 55, 333–361.
32. J. Kästner, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 932–942.
33. F. Schmitz, P. Virnau and K. Binder, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 053302.
34. L. Dagum and R. Menon, IEEE Comput. Sci. Eng., 1998, 5, 46–55.
35. R. J. Allen, C. Valeriani and P. R. ten Wolde, J. Phys.: Condens. Matter, 2009, 21, 463102.
36. F. A. Escobedo, E. E. Borrero and J. C. Araque, J. Phys.: Condens. Matter, 2009, 21, 333101.
37. R. J. Allen, P. B. Warren and P. R. ten Wolde, Phys. Rev. Lett., 2005, 94, 018104.
38. S. Ružička, D. Quigley and M. P. Allen, Phys. Chem. Chem. Phys., 2012, 14, 6044–6053.
39. J. S. Langer, Phys. Rev. Lett., 1968, 21, 973–976.
40. V. A. Shneidman, K. A. Jackson and K. M. Beatty, J. Chem. Phys., 1999, 111, 6932–6941.
41. G. Jacucci, A. Perini and G. Martin, J. Phys. A: Math. Gen., 1983, 16, 369–383.

---