Impact of Surface Roughness on Crystal Nucleation

Abstract

We examine the effect of rough surfaces on crystal nucleation by means of kinetic Monte Carlo simulations. Our work makes use of three-dimensional kMC models, explicit representation of transport in solution and rough surfaces modeled as randomly varying height fluctuations (roughness) with exponentially decaying correlation length (topology). We use Forward-Flux Sampling to determine the nucleation rate for crystallization for surfaces of different roughness and topology and show that the effect on crystallization is a complex interplay between the two. For surfaces with low roughness, small clusters form on the surface but as clusters become larger they are increasingly likely to be found in the bulk solution while rougher surfaces eventually favor heterogeneous nucleation on the surface. In both cases, the rough surface raises the local supersaturation in the solution thus leading to another mechanism of enhanced nucleation rate.

1. Introduction

Crystallization is of broad importance in many areas of science and plays a key role in numerous industrial processes such as the production of pharmaceuticals. However, many molecules of interest are difficult to crystallize, particularly macromolecules such as proteins and monoclonal antibodies. This has long motivated the search for means to enhance crystallization via the use of seeds [1] and substrates, or surfaces, References [2,3] that can induce heterogeneous nucleation of the crystalline phase. Important theoretical insight was gained from the seminal study of Page and Sear [4] showing that porous substrates—with pores many times the size of a molecule—enhance crystallization more than flat substrates (see also [5]). This work was extended to rough substrates in which the critical nucleus is much larger than the length scale of the substrate inhomogeneity [6]. Experimental work has also demonstrated the value of rough substrates as nucleants [7,8] while the importance of adsorption of biomolecules such as proteins on rough substrates is important in applications such as bioceramics [9] and in controlling the stability and safety of proteins in pharmaceutical formulations [10]. Despite such motivation, existing theoretical studies have severe limitations that make their translation to real systems unclear. Our goal here is to revisit the problem of rough substrates using similar, but more realistic, techniques.

In the simulation studies mentioned above, standard two-dimensional Ising models were used in which space is divided into a square lattice and each lattice site possesses a label or “spin” that can take on two values: e.g., “liquid” or “solid”. Two adjacent cells with same label contribute an energy $-J<0$ to the system energy while two with opposite labels contribute $+J>0$ to the energy. Supersaturation is introduced by a contribution of $-h<0$ for each cell labeled with the nucleating phase (say, the solid phase) thus making it preferred. Walls are cells containing no spins and not contributing to the energy of any adjacent cells. Initially, all of the cells are in the mother phase (e.g., “liquid”) and the system is evolved by selecting on cell at random and changing its label according to the standard Metropolis algorithm [11]. Rare event techniques, such as Forward Flux Sampling (FFS) [12,13], were used to simulate nucleation and it was found that nucleation occurred much faster in the presence of pores or rough substrates.

These types of Ising models lack several key ingredients of the physics of real systems. Most importantly, many of the applications such as protein crystallization involve nucleation from a low-concentration (or weak) solution. This means that the formation of a dense solid nucleus requires the molecules in solution over a large volume to come together and this mass transport is, therefore, a limiting factor in the dynamics and thus plays a central role in the process [14,15]. Since the standard Ising model has nothing like the conservation of mass, this important physical element is completely absent. Once this is added to the model in the form of molecules that move around, dimensionality also becomes important so that one must be cautious in extrapolating 2D simulations to the real 3D world.

Here, we try to improve upon these points in the obvious way: we use three-dimensional models with mass conservation and a realistic diffusional dynamics in the solution. As in the previous studies, we use FFS to allow us to track nucleation. Furthermore, unlike the work of [3,6] we use standard models of rough surfaces that realistically model both the height probability distribution of the surface and its autocorrelation. We find that molecular adsorption of rough substrates to be enhanced compared to flat substrates which can be interpreted as arising due to the greater surface area of a rough surface as well as the fact that molecules can interact with more than one surface site at a time. However larger surface area alone does not explain the enhanced molecule/substrate interaction. Our results show that increasing surface roughness leads to localization of the molecules on the substrate. We also find that surface correlation length plays an important role in determining system properties like free energy (enthalpic effect) and in the adsorption/desorption transition. Of most relevance here, we observe a transition in the nature of crystal growth for critical values of the roughness and correlation length. In summary, our results show that substrate roughness enhances crystallization rates but that the process is more complex than simply heterogeneous nucleation on the substrate. We find that the dominant effect is the enhancement of the local concentration of solute near the substrate and that it is via this indirect mechanism that crystal nucleation rates are enhanced.

In the next section, we present our simulation model: the model for the solution, which is essentially a standard Ising-type kinetic Monte Carlo model with diffusional dynamics, the modeling of surface roughness and the determination of nucleation rates. In the following section, we present results and analysis focusing on the nucleation rate for substrates with different roughnesses. In particular, we try to unravel the dual effects of roughness in terms of the amplitude of height variations and of the correlation of height variations over the substrate. Finally, we close with a discussion of the implication of our results for controlling crystallization in real systems.

2. Computational Methods

2.1. Simulation Model

In our model, the simulation cell is partitioned into a lattice of $N=N_x\times N_y\times N_z$ cubic cells with sides of length a, the typical size of a molecule. The lower part of the lattice is occupied by the substrate: a random surface described mathematically as $h=h(x,y)$, where h is the surface height with respect to bottom ($z=0$) of the lattice, and $(x,y)$ is the position vector on the surface. Above the substrate, the lattice sites are either empty or occupied by a single molecule. The total energy is determined by counting the number of nearest-neighbor bonds with each molecule-molecule bond contributing energy $-ϵ$ and each molecule-wall bond contributing energy $-{ϵ}_W$ where, in our simulations, both $ϵ,{ϵ}_W>0$.

Let the total number of mobile molecules (i.e., those with one or more unoccupied nearest neighbor site) in the simulation a given time t be $N_t$. Then, a kMC move consists of randomly choosing a number, call it n, between 0 and $N_t+N_x\times N_y$ and attempting to execute one of two possible actions. If $n<N_t$ then an attempt is made to move molecule n to one of the neighboring lattice sites: a site is chosen at random and the move is accepted with probability $min(1,e^{-\Delta E/k_{B}T})$ where $\Delta E$ is the change in energy of the system if the move is accepted. In the x- and y-directions we apply periodic boundary conditions while in the z-direction molecules can exit the simulation cell. These nearest-neighbor hops simulate molecules diffusing in the solution. If the number chose is $n>N\left(t\right)$ then an attempt is made to add a molecule to the cell $n-N_t$ in the top layer of the system—the probability that a molecule appears is determined by an input parameter. In either case, the time is advanced by a variable time step, $\Delta t_t={\nu }_0/N_t$ where the “attempt frequency” ${\nu }_0$ is a (fixed input) parameter having units of inverse time. The physical value of the attempt frequency can be determined by considering the tracer diffusion of a single molecule in solution. For a single molecule, the kMC algorithm generates a discrete random walk with diffusion constant $D=a^2{\nu }_0$ so that in principle, real systems can be mapped to the model by taking the lattice spacing a to be the typical size of a molecule and then choosing the attempt frequency ${\nu }_0$ so as to give the observed single-molecule diffusion constant. However, for simulations, one takes a to be the length scale and ${\nu }_0$ the time scale and so the physical values are not needed.

Since molecules can leave the simulation cell by hopping out at the top and are periodically inserted, an equilibrium is reached. Inside the simulation cell, mass is conserved by the diffusive dynamics, whereas the upper boundary represents an interaction with a reservoir. This simulates an open system and is necessary to simulate nucleation of a dense phase from a weak solution: otherwise, the material in a closed simulation cell would be depleted in forming the cluster.

2.2. Model of Rough Surface

Many roughness effects are closely related to the interplay between the length scales characterizing the roughness and the adsorbing molecule [16]. In the case of macromolecules, for example, due to their complex structure, entropic as well as energetic effects play a role in determining their binding to a surface [17]. Length scale matching is particularly important for interactions involving hydrophobic parts of the surface of the protein since roughness can change surface hydrophobicity. Since in our model, a molecule occupies a single lattice site, the lattice spacing a is taken to be the typical size of a molecule and we note that molecules making up the substrate and those of the crystallizing material are therefore of the same size.

In general, a surface is fully specified by giving its height $h\left(\mathbf{r}\right)=h(x,y)$ for each pair of coordinates $x,y$ or, for our discrete model, one could also write $h(x,y)=h(ka,la)\equiv h_{kl}$ for some integers $k,l$. A rough surface corresponds to one for which the heights $h_{kl}$ are random variables with some specified statistics. The variability of surface height is characterized by the probability $p\left(h\right)dh$ that the height at any given cell is between h and $h+dh$, i.e., the height distribution function for which we use a standard Gaussian model with variance $\sigma$,

$$p\left(h\right)=\frac{1}{\sqrt{2\pi }a\sigma }{e}^{-\frac{h^2}{2a^2{\sigma }^2}}.$$

(1)

This does not, however, fully characterize the surface since there can be height correlations so that two cells that are close together have a higher (or lower) probability of having similar heights than two which are far apart, which is described by the height correlation function $G({\mathbf{r}}_1,{\mathbf{r}}_2)=\langle h\left({\mathbf{r}}_1\right)h\left({\mathbf{r}}_2\right)\rangle$. In principle, the average in this expression should be understood as an ensemble average: an average over many different realizations of the surface. In the following, we restrict attention to homogeneous, isotropic rough surfaces for which the height correlation function depends only on the distance between positions ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$, e.g., $G({\mathbf{r}}_1,{\mathbf{r}}_2)=G\left(\right|{\mathbf{r}}_1-{\mathbf{r}}_2\left|\right)=G\left(r_{12}\right)$. In this case, a measure of the typical length scale of such correlations, i.e., of the lateral correlation length, is given by ${\int }_0^{\infty }G\left(r\right)rdr/{\int }_0^{\infty }G\left(r\right)dr$. In our work, we have used the standard exponential model,

$$G\left(r\right)=a^2{\sigma }^{2}exp(-r/\left(a\xi \right)),$$

(2)

where it can be confirmed that the parameter $\xi$ corresponds to the correlation length (divided by the lattice spacing) defined above. Our goal is to generate surfaces—realizations of the stochastic process—for which, in the infinite-surface size limit, the heights reproduce the chosen height distribution and correlation function. To do this, we use the Fourier transform method described in Ref [18], details of which are given in the Appendix A.

Following Wenzel [19,20], we introduce a more physical measure of the roughness of the surface $r_w=A/A_0$, known as the Wenzel roughness or Wenzel factor, $r_w$, which is the ratio of the surface area A of the rough surface ot its projected area $A_0$. For the exponential model, it is shown in the Appendix A that this is given by the sum

$$\overline{r_w}\simeq 1+2\frac{{\sigma }^2}{\xi }-4\frac{{\sigma }^4}{{\xi }^2}+8\frac{{\sigma }^6}{{\xi }^3}+\cdots$$

(3)

and in the high roughness limit, one gets

$$\overline{r_w}=1+\frac{4}{\sqrt{\pi }}\sigma \sqrt{1-e^{-1/\xi }}.$$

(4)

2.3. Rate Calculations

In this study, we compute heterogeneous nucleation rates using the Forward Flux Sampling algorithm (FFS). We identify a cluster as being a set of molecules connected by nearest-neighbor bonds and we use the size of the largest cluster at any given moment, $\lambda$, as the order parameter for the system. In the initial, metastable, state consisting of a weak solution of crystal-forming molecules, only small clusters are commonly observed. We, therefore. characterize it as having $\lambda \le {\lambda }_{\mathrm{S}}$ where ${\lambda }_S$ is a somewhat arbitrary threshold. The post-critical, crystalline state is characterized by having a sufficiently large cluster (i.e., one much larger than the critical cluster) and so can be defined as ${\lambda }_X\ge {\lambda }_{\mathrm{C}}$ for some (crystal) threshold ${\lambda }_{\mathrm{C}}$.

The nucleation rate coefficient $k_{SC}$, considered as the frequency of forming one stable large cluster from the random solution, is given by the cumulative flux of trajectories that begin with $\lambda <{\lambda }_{\mathrm{S}}$ and end up in $\lambda >{\lambda }_{\mathrm{C}}$ and is estimated as the product of ${\Phi }_0$, the initial flux of trajectories that cross ${\lambda }_{\mathrm{S}}$ and of the probability that a system that has $\lambda ={\lambda }_{\mathrm{S}}$ evolves into one with $\lambda ={\lambda }_{\mathrm{C}}$ without ever returning to the initial state. To determine this probability, a system in the solution state is allowed to evolve and each time, it crosses from $\lambda <{\lambda }_{\mathrm{S}}$ to $\lambda ={\lambda }_{\mathrm{S}}$, the configuration is stored and the simulation continues. In this way, a database of crossings is collected and, at the same time, the number of crossings per unit time—i.e., the flux—is determined. Next, an intermediate threshold ${\lambda }_1>{\lambda }_S$ is defined and the probability of the transition ${\lambda }_S\to {\lambda }_1$ is determined by randomly choosing a configuration from the database of configurations with $\lambda ={\lambda }_S$ and evolving it until it either reaches ${\lambda }_1$, in which case it is stored in a new database, or falls back into the metastable state, in which case it is discarded. The fraction reaching ${\lambda }_1$ gives the probability of the transition. This is then repeated for a series of such milestones, ${\lambda }_{\mathrm{S}}<{\lambda }_1<{\lambda }_2<\cdots <{\lambda }_N={\lambda }_{\mathrm{C}}$ and the various transition probabilities are determined. The product of all of them then gives the overall transition probability and, when multiplied by the flux out of the initial state, gives the transition rate and this divided by the volume gives the nucleation rate. In our simulations, we took ${\lambda }_{\mathrm{S}}=5$, ${\lambda }_{\mathrm{C}}=100$ and milestones at each integer in between. From each milestone ${\lambda }_k$, new trajectories are launched until there are 400 that reach the next milestone.

3. Results and Discussion

We performed kinetic Monte Carlo simulations of the lattice model on a system with size $N_x=N_y=N_z=50$ in lattice units a. In all the simulations the temperature is fixed at $k_{B}T=0.4ϵ$ and the bulk density fixed by the reservoir is $\rho =0.0016$ in lattice units. The bulk density was chosen sufficiently low that the nucleation was not instantaneous—making impossible to meaningfully study the effect of substrates—but sufficiently high that the critical cluster was not large compared to system dimensions. (Typical critical clusters consisted of between 60 and 100 molecules, so with diameters between 5 and 6 lattice sites.) We generated exponentially correlated random rough surfaces with given variance ${\sigma }^2$ and correlation length $\xi$, typical examples of which are shown in Figure 1.

3.1. For Hard Walls, Roughness Is not Important

The first question addressed was whether purely entropic—i.e., geometric—effects, induced by the substrate topography may influence the nucleation rate through an effect of confinement.

Figure 2 shows the rate of formation of clusters when the solution is in contact with purely repulsive substrates (${ϵ}_W=0$) or substrates slightly attractive (${ϵ}_W=0.1$) with varying topographies. As can be seen, the rate of formation of clusters is unaffected by changes in the variance of the height indicating that entropic effects alone, in the absence of significant attractive interactions, do not significantly affect the nucleation rate.

In contrast, when the substrate is attractive (${ϵ}_W=0.3$), the nucleation rate becomes quite sensitive to the surface topography. The combined effect of enthalpy and surface roughness is clearly seen in Figure 3 showing the critical sizes of the clusters and Figure 4 which shows the rate of cluster formation for different values of the height variance with correlation lengths $\xi =5$ (upper graphic), and $\xi =20$ (lower graphic). For the different surface topographies investigated, at fixed correlation length, the critical size decreases from about 87 to just above 60 and the nucleation rate increases by a factor ${10}^7$ as the variance increases from ${\sigma }^2=0.5$ to ${\sigma }^2=4.5$. Moreover, at fixed variance the nucleation rate increases by a factor of a thousand as the correlation length decreases from $\xi =20$ to $\xi =5$. The critical size decreases almost monotonically for the smaller correlation length but at larger correlation length it is unaffected until the roughness exceeds ${\sigma }^2=4$. Our goal in the following is to understand the mechanisms behind these variations.

3.2. Roughness of Attractive Walls Increases Absorption

The usual explanation for the variation of nucleation rate with roughness is that a rough substrate has a larger surface area than a flat substrate so that more molecules can adsorb onto it leading to heterogeneous nucleation. To understand the effect of roughness on molecular adhesion one can follow the thermodynamic analysis and assume that the amount of molecules in contact with the substrate is proportional to the energy of adhesion $w_{ad}$ which for a flat solid surface is defined by [21,22]

$$w_{ad}=\left({\gamma }_{cv}+{\gamma }_{sc}\right)-{\gamma }_{sv},$$

(5)

where ${\gamma }_{sv}$, ${\gamma }_{sc}$, and ${\gamma }_{cv}$, are the interfacial tension for substrate–vapor, substrate–crystal, and crystal–vapor interfaces. By combining (5) with Young’s equation ${\gamma }_{sv}-{\gamma }_{sc}={\gamma }_{cv}cos{\theta }_0$ one obtains the Young–Dupre equation

$$w_{ad}={\gamma }_{cv}\left(1-cos{\theta }_0\right),$$

(6)

where ${\theta }_0$ is the equilibrium contact angle of the crystal on the substrate. For a rough surface the adhesion energy and the contact angle become spatial-dependent [23]

$$w_{ad}\left(\mathbf{r}\right)={\gamma }_{cv}(1-cos\theta \left(\mathbf{r}\right)),$$

(7)

where $\theta \left(\mathbf{r}\right)$ is the contact angle modified by the local slope at the surface, $\theta \left(\mathbf{r}\right)={\theta }_0+arctan\left(\right|\nabla h\left|\right)$. In the discrete model, there is a perfect contact between molecules and the surface since the lattice spacing is identical to the molecular size. In this case the average energy of adhesion is proportional to number of plaquettes of the solid surface in contact with molecules, $w_{ad}={ϵ}_W〈A_d〉$, and using Equation (A25) from the Appendix A, gives the high-roughness estimate

$$w_{ad}={ϵ}_W\left(1+\frac{4}{\sqrt{\pi }}\sigma \sqrt{1-e^{-\frac{1}{\xi }}}\right).$$

(8)

This shows that for fixed correlation length, adhesion energy increases with increasing variance, $\sigma$. Conversely, for fixed variance, adhesion energy decreases with increasing correlation length. All of this is consistent with the intuition that crystal that forms in small valleys on the substrate will have more contact with the substrate and so lower energy than the same volume of crystal on a flat substrate: in other words, roughness over molecular length scales lowers the crystal–substrate surface tension. Assuming that the amount of absorbed molecules at the substrate is proportional to the adhesion energy, the effect of roughness on adhesion is demonstrated by measuring the molecular surface density ${\rho }_s$ as a function of $\sigma$ and $\xi$. The surface density ${\rho }_s$ is defined as the number of molecules whose distance from the plane $z=0$ is lower so that $h_{\max }=max\left(\{h\left(\mathbf{r}\right):\mathbf{r}=(x,y)\}\right)$ is normalized by the surface area of the plane, $N_x\times N_y$.

Figure 5 shows that ${\rho }_s$ is indeed proportional to $\sigma$, and that by increasing the correlation length one reduces the surface density. Furthermore, the observed local increase in concentration is well-correlated with the nucleation rate as shown in Figure 6.

3.3. Correlation Plays an Important Role

While these results ascertain that rougher substrates adsorb more molecules and that this is correlated with a higher nucleation rate, Figure 5 suggests that ${\sigma }^2$ is not the only factor. Clearly, substrates with the same roughness but different correlation lengths absorb different numbers of molecules so that the topography also plays an important role.

A clue as to the role of surface correlations comes from an examination of post-critical clusters, see Figure 7. Contrary to what one might expect, clusters are not all in contact with the substrate and, in fact, some grow far away from it. To quantify this, we measured the contact surface $A_{sc}$, between the substrate and the largest cluster, corresponding to $\lambda =100$, as a function of ${\sigma }^2$ and $\xi$.

The probability distribution of $A_{sc}$, obtained from a thousand substrate–cluster configurations, is displayed Figure 8. One observes that for ${\sigma }^2=2.5$ the most probable configuration is that the cluster has no contact with the substrate while at higher variances, the probability of no-contact decreases until at ${\sigma }^2=5$, all clusters grow on the solid substrate. To understand this better, Figure 9 shows the evolution of the contact area for different values of ${\sigma }^2$, as a function of the size of the growing crystal: in other words, how the contact area changes as a function of cluster size.

The behavior of $P\left(A_{sc}\right)$ as a function of cluster size demonstrates that when the roughness is small, crystals start to grow on the substrate that favors heterogeneous nucleation, but they detach from the substrate when reaching a critical size ${\lambda }_c$, while for rougher substrates, they always remain attached to the substrate. A similar behavior is observed in seeded heterogeneous crystallization [24]. In [24], the sequence of crystallite structures is followed during crystallization in a colloidal system in contact with a spherical glass bead. Depending on the bead diameter, the crystallite detaches from the curved seed surface once it reaches a critical size. The origin of this behavior is attributed to a mismatch between the crystalline structure and the curved surface that induces elastic stress in the crystal. Because the distortion that leads to the detachment depends on the crystallite-seed relative curvature, the critical crystal diameter is a function of the seed diameter. In the present model, the detachment could not be attributed to the relaxation of an elastic distortion, however, a similar conclusion can be drawn by considering the interfacial energy of the crystal.

To begin with, it is clear that a cluster of two molecules will be more stable if it is in contact with the substrate (which contributes some bonding energy) than in the bulk fluid. For this reason, small clusters are more likely to be found in the valleys of the rough surface where they can maximize their contact with the substrate (or, in other words, replace as much vapor–crystal surface with crystal–substrate surface as possible). For larger clusters, there are two competing effects. First, if they form in the valleys of the substrate, they continue to replace the crystal–vapor interface with the crystal–substrate interface, which is energetically favored. Indeed, one can estimate the crystal–vapor surface tension to be $\left|ϵ\right|/2a^2$ whereas the crystal–substrate surface tension is ${\gamma }_{\mathrm{CS}}=|(ϵ/2)-{ϵ}_W|/a^2$. Second, however, is that in bulk, a cluster can minimize its surface area (by becoming spherical) and hence lower its energy: this is not possible for clusters growing in the substrate which are constrained by the geometry of the pores. Thus, one expects that small clusters will be energetically favored to form on the substrate while large clusters may have lower energy in the bulk, away from the substrate. An important question here will be the stability of the clusters that form in the pores. Small clusters will, generally, be unstable but there may be a critical size beyond which the clusters in the pores become stable. In this case, one expects to find large clusters growing on the substrate. However, as the clusters grow they fill the pores and it may be that stability is never attained until they exceed the size of a pore in which case, they may be less stable than equal-sized clusters in the bulk.

This picture is confirmed by Figure 9 which shows the average area of contact with the substrate as a function of cluster size. Note that this is an average over all clusters including those that are not in contact with the substrate at all. As expected from our argument, for the smallest clusters, the average contact area increases with increasing cluster size regardless of the roughness. For the case of low variability of the surface height, a maximum is reached at a certain cluster size and the contact area then decreases until eventually reaching zero, indicating that detached clusters become increasingly favored. For the rougher cases, however, the contact area increases monotonically showing that clusters remain attached to the substrate as they grow. Figure 10 shows the average volume of a valley on the substrate as a function of variance of the height. We find that, as one would expect, the volume of the valleys increases with roughness and, furthermore, that the volumes are close to the observed maxima of the contact area at lower roughness. This suggests that bulk clusters become more stable than surface clusters once the clusters fill the pores unless, as argued above, the pores exceed a critical volume in which case clusters on the substrate can continue to grow indefinitely.

The overall impact of the substrate roughness on the nucleation rate depends on both the surface geometry (i.e., its roughness is controlled by both the variance of the height and the correlation length) and on the substrate chemistry (here characterized by the binding energy of the crystal molecules to the substrate). Clearly, as the strength of the crystal–substrate bond decreases, the crystal–substrate surface tension increases and the geometric constraint of the pores makes it increasingly unlikely that clusters will form in them. This is illustrated in Figure 11 which shows the nucleation rate as a function of crystal–substrate bond strength for a fixed variance (${\sigma }^2=2.5$) and two different correlation lengths. In both cases, the nucleation rate is unaffected by the roughness up to a critical bond strength beyond which the substrate enhances the nucleation rate. As expected, the enhancement is greater for the smaller correlation length since in that case, one expects the pores to have a larger average volume.

This scenario depends on both the length scales, and its impact on the nucleation rate, as explained above, is heavily dependent on the binding energy. Figure 11 shows the nucleation rate as a function of the binding energy for two values of the correlation length. When ${ϵ}_W<{ϵ}_c$, where ${ϵ}_c$ is a critical binding strength, the nucleation rate is almost independent of ${ϵ}_W$, while for ${ϵ}_W>{ϵ}_c$, the nucleation rate increases with ${ϵ}_W$. However, increasing interaction leads to irreversible adsorption on the substrate, therefore hence the limit $|{ϵ}_W|\le 0.3$.

4. Conclusions

Our work reveals fundamental insights into the role of roughness and chemistry on nucleation and, in particular, on the complex interplay of the variance of height of the surface and the correlation length. We identify two important physical phenomena that contribute to the observed substrate-enhanced nucleation. The first is that any substrate that forms bonds with the crystal molecules enhances the local concentration of the crystal molecules near the surface and this correlates with increasing nucleation rate, regardless of details of the substrate geometry. In this effect, the main role of roughness is that it increases the effective surface area and so increases the substrate-generated supersaturation with respect to the bulk fluid.

A second effect is that small clusters preferentially form in the valleys of the rough surface. When the crystal–substrate bonding is sufficiently strong, cluster formation within the pores is favored over that in the bulk at least until the clusters fill the pores. If the pores are small (i.e., the roughness is low) then larger clusters form preferentially in the bulk but if the pores exceed a critical size, then the growth of the pore-nucleated clusters continues and is the dominant mechanism for the formation of bulk crystal.

Our work has a number of practical, engineering-related implications. While one usually imagines controlling surface absorption by changing the substrate chemistry, our results indicate that even for fixed substrate chemistry, the surface absorption—and, hence, nucleation rate—can be tuned by changing the topography in the form of the correlation length. Our results indicate that less correlated surfaces have higher surface absorption. This could enable substrate engineering in situations where chemical modification is not possible. Furthermore, we have shown a clear mechanism for understanding substrate-dominated crystallization versus substrate-enhanced crystallization in the bulk which could prove useful, e.g., in engineering membranes for membrane-enhanced crystallization. Finally, we note that these observations are only possible in realistic models that include the bulk fluid, concentration gradients and realistic models of roughness. Studies based on simple solid-on-solid substrate models do not include any of these elements and so can not fully reveal the complexity of the process.