Mechanism of two-step vapour–crystal nucleation in a pore

We present a numerical study of the effect of hemispherical pores on the nucleation of Lennard–Jones crystals from the vapour phase. As predicted by Page and Sear, there is a narrow range of pore radii, where vapour–liquid nucleation can become a two-step process. A similar observation was made for different pore geometries by Giacomello et al. We find that the maximum nucleation rate depends on both the size and the adsorption strength of the pore. Moreover, a poe can be more effective than a planar wall with the same strength of attraction. Pore-induced vapour–liquid nucleation turns out to be the rate-limiting step for crystal nucleation. This implies that crystal nucleation can be enhanced by a judicious choice of the wetting properties of a microporous nucleating agent.


Introduction
Freezing in three-dimensional systems is a first-order phase transition. As a consequence, it is possible to supercool or super-compress the liquid parent phase until crystallisation happens spontaneously on the timescale of the experiment. Often the crystallisation process starts with a nucleation step during which a crystalline nucleus is formed that can then grow to macroscopic size, and the overall rate of crystallisation is determined frequently by the nucleation step. Moreover, the structure of the crystal that grows to macroscopic size is often determined by the structure of the crystal nucleus that forms during the rate-limiting step. For this reason, it is of both fundamental and practical interest to gain insight into the mechanisms that can be used to control the rate and the pathway of crystal nucleation.
Nucleation can proceed either homogeneously, that is, in the bulk of the parent phase, far away from walls or other external influences, or heterogeneously at the surface of a 'seed'. In most cases of practical interest, heterogeneous nucleation dominates the rate of crystal formation [1].
For this reason, heterogeneous crystal nucleation has been studied extensively ever since the early work of Ostwald [2]. In recent years, there has been something of a revival of studies of heterogeneous nucleation both from a theoretical/numerical perspective and from experiments. This renewed interest was stimulated, on the one hand, by the availability of suitable colloidal model systems that make it possible to follow the dynamics of individual particles during nucleation in real space and real time [3,4] and, * Corresponding author. Email: df246@cam.ac.uk † Current address: Connex.cc DI Hadek GmbH, Nordbahnstrae 56/DG2, A-1020 Vienna, Austria. on the other hand, by the development of novel computational techniques that made it possible to study the early stages of crystal nucleation in great detail [5,6].
Experiments by Chayen and co-workers [34,35] on the crystallisation behaviour of a number of different proteins demonstrated that the rate of protein crystallisation can be strongly enhanced in the presence of a microporous medium. This intriguing observation stimulated subsequent theoretical work by Page and Sear [25], who studied the nucleation behaviour of a (two-dimensional) Ising model in a rectangular indentation in a wall. This lattice model did indeed show that pores can enhance nucleation and that this effect depends strongly on the pore diameter.
The key findings of Page and Sear [25] were reproduced in the subsequent lattice density functional theory of Liu et al. [27]. More recently, extensive numerical simulations of Giacamello et al. [28] have further elucidated the C 2015 The Author(s). Published by Taylor & Francis. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/Licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. mechanism of two-step vapour-liquid nucleation for a variety of pore geometries. Turning now to crystal nucleation, it is clear that a lattice model such as the one used in [25] cannot capture one of the most interesting features of protein crystallisation, namely that it is often most effective in a temperature-concentration range where the protein solution is close to liquid-liquid demixing [36][37][38]: the Ising model does not distinguish between the liquid and the crystalline state. It is, therefore, of interest to expand the pore nucleation studies to an off-lattice model for which the liquid and the crystalline states are distinct. With such a model study, we can look in more detail at the ways in which a porous medium affects the pathway for crystallisation. In particular, we can address questions such as: does a crystal form immediately in the pore or does crystallisation proceed via an intermediate liquid phase? Another advantage of an off-lattice model is the fact that the diameter of the pore needs not be commensurable with the lattice spacing of the crystal (in a normal lattice model, this commensurability is inevitable). If a crystal does not 'fit' optimally in a pore, one might expect that this could affect the free-energy barrier for crystal nucleation.
In an earlier publication [39], we provided evidence that capillary condensation in a pore may speed up crystal nucleation. The present paper describes a systematic study of the effect of pore shape, size and interaction strength.
In the present study, we first investigate the effect of a hemispherical pore on the nucleation pathway for vapourliquid nucleation. In Section 3.1, we explore nucleation at a planar (weakly) attractive wall. Next, we consider nucleation on a planar circular attractive patch, and finally, we consider a pore. Section 3.2 focuses on liquid-crystal nucleation. We find that, as in the case of homogeneous nucleation, this is an independent nucleation event that can be treated separately. Finally, Section 4 integrates the resulting picture for the overall nucleation pathway in systems with surface heterogeneities.

Simulation details
The present study is based on Metropolis Monte Carlo simulations of a model system consisting of particles interacting through a truncated and shifted Lennard-Jones (LJ) pair potential [40]: where the full (i.e. not truncated) LJ interaction is given by Here is the unit of energy, σ is the unit of length and r c the interaction cut-off distance. Note that the choice of interac-r c R Figure 1. Schematic side view of a pore. The dark grey area indicates the space that is inaccessible to LJ particles due to the presence of a hard wall. The range of the weak attraction is depicted in light grey. Molecules in the white area are outside the range of interaction of the wall. The radius of the pore is denoted by R and the range of attraction of the pore walls is denoted by r c .
tion cut-off has a significant effect on the free energy of the system and on properties such as coexistence curves and surface tensions [41,42]. However, the qualitative features of the phase diagram remain unaffected. Hence, the model that we use is adequate for the present study that focuses on generic effects, rather than on the properties of a specific experimental system. Three different geometries were used to study heterogeneous nucleation: a planar wall, a planar circular patch and a hemispherical pore. Both the circular patch and the hemispherical pore are embedded in an otherwise purely repulsive substrate. Figure 1 depicts a schematic side view of such a pore, which consists of an open half-sphere embedded in a hard wall. It shows the relevant parameters: the pore radius R and the pore's interaction range r c . All surface geometries are weakly adsorbing due to an effective potential given by the LJ 9-3 potential, which results from integration of the pair interactions over the half-space, assuming a constant particle density. Here, r is defined as the distance perpendicular to the surface. The potential U w (r) was cut and shifted at a distance r c , and we choose r c and σ to be the same for both U w (r) and U(r).
In what follows, we use reduced units. We define the reduced distance as r * = r/σ and the reduced potential energy as u * = U/ . All other reduced quantities (e.g. the pressure P * = Pσ 3 −1 , the density ρ * = ρσ 3 and the temperature T * = k B T −1 ) follow. All quantities reported in this work are expressed in these reduced units. Therefore, we omit the superscript asterisk (*) from here on.
We performed all vapour-liquid simulations in the grand canonical ensemble, where temperature, volume and the chemical potential are kept constant, and the particle number is allowed to fluctuate. This ensures a constant vapour pressure and minimises finite size effects. In the simulations of liquid-solid nucleation, we consider a liquid droplet in contact with both the vapour phase and the surface. These simulations were performed in the canonical ensemble to keep the size of the liquid droplet constant. The initial liquid droplet was taken from grand canonical vapour-liquid simulations, then equilibrated at a high temperature to remove all crystallinity, whilst constraining the system (using umbrella sampling) to prevent evaporation. The droplet was then quenched to T = 0.45 and equilibrated using umbrella sampling to keep the droplet in the metastable liquid phase. This procedure allowed us to produce metastable liquid droplets of the required size.
Periodic boundary conditions were applied in the y-and z-directions of the cubic simulation boxes, and a planar hard wall prevented particles from escaping in the x-direction. For all simulations, we chose the box size to be l 3 = 50 3 and the particle displacement trial move size to be x = 0.2. In the grand canonical ensemble, each particle was moved on average 20 times before an insertion/removal move was performed.

Rare event sampling
In order to study the properties of the system on the path towards nucleation, we applied an umbrella sampling [43] (US) scheme, where a biasing potential is used to constrain the size of the nucleus to be close to a certain size n 0 . By computing the bias-corrected cluster size probability distribution P 0 (n) for many different n 0 , a series of overlapping windows is obtained. By fitting these windows together, the free energy as function of cluster size can then be computed from the overall probability histogram, G(n) = −k B T ln P(n) + c, where c is a constant offset. To be more precise, we define P(n) as the ratio of M(n), the number of clusters of size n, divided by the (average) total number of particles N. For large cluster sizes, P(n) 1 and the probability of observing more than one cluster of this size (or of a larger size) in a simulation is negligible. In that case, the probability to find a cluster of size n in the system can be approximated by the probability to find that the largest cluster in the system contains n particles. This facilitates the umbrella sampling because we can then apply a bias that controls only the size of the largest cluster in the system [6]. However, when computing P(n) for small clusters, we cannot use this approach as there will be many clusters with sizes comparable to, or larger than n. In that case, we sample P(n) from an unbiased run taking all cluster sizes into account. This is important because if we would constrain the size of the largest cluster for small n (thereby cutting off the cluster distribution for all sizes larger than n), we would be truncating the natural cluster size distribution in the system. As a consequence, the free energy of the system would then appear to increase as n → 1, and as a result, there would be a spurious (and system-size-dependent) minimum in the free energy G(n) at small n.
For the computation of the vapour-liquid nucleation free-energy barriers, we started from the pure vapour phase. For the liquid-crystal nucleation, we first equilibrated a liquid droplet of n l ≈ 700 particles, biasing the system to prevent any crystallisation. Then, we performed US simulations in the constant number, volume and temperature (NVT) ensemble to compute the nucleation free-energy barrier.

Order parameters
For the vapour-liquid nucleation, we used an order parameter based on Stillinger's overlapping spheres criterion [44]. It defines a particle to be in a high-density phase, if it has at least one neighbour within a distance 1.5 corresponding roughly to the first minimum of the pair correlation function of a bulk liquid. Then, a cluster analysis is performed on all high-density particles and the number of particles in the largest cluster is taken as the order parameter. Note that this order parameter does not distinguish between an ordered or disordered high-density phase, and therefore, does not favour one phase (liquid or crystal) over the other. For a detailed description of this order parameter, we refer the reader to [45,46].
As order parameter for the liquid-solid nucleation, we applied the local bond-order parameter [6,47]. This order parameter assigns each particle a 13-dimensional vector capturing its local environment: where Y 6m denotes a sixth-order spherical harmonic with components m ranging from −6 ≤ m ≤ 6. N b (i) is the number of nearest neighbours of particle i, andˆ r ij a unit vector connecting the centres of mass of particles i and j. The sum is overall neighbouring particles j within a cut-off distance r = 1.5. In the second step, the order parameter computes the dot product q 6 (i) · q 6 (j ) between each particle i and all its neighbours j, effectively comparing the particles' neighbourhoods. If q 6 (i) · q 6 (j ) exceeds a threshold of 0.65, particles i and j are considered to form a 'link'.
Only if a particle's total number of links n(i) exceeds five, it is considered to be solid-like. In the final step, Stillinger's criterion is applied to all solid-like particles to identify the number of particles in the largest solid-like cluster, which is used as order parameter. All parameters for this order parameter can be obtained by minimising the overlap between distributions of an equilibrated bulk solid and a metastable bulk liquid at working conditions. We note that neither the nucleation rate nor the nucleation pathway of our simulations was sensitive to (moderate) variations in the definition of the order parameters. However, such changes do affect the apparent size of the clusters that we detect, including the critical cluster size. The local bond-order parameter that we used was originally designed to detect nucleation in the bulk. In order to verify that nucleation at the droplet surface is still properly detected, we also performed simulations with the modified local bond-order parameter used by Mendez-Villuendas and Bowles [48]. This modification ensures that surface particles, too, can be identified as solid-like particles and are taken into account properly. However, no significant differences were observed.
For all simulations involving a heterogeneity, the order parameter is adapted to identify only clusters in contact with the heterogeneity. This means that in addition to the procedure mentioned above, only those clusters are considered that contain at least one particle that is within the effective interaction range of the surface heterogeneity. We refer to this modification as a restricted order parameter (rOP). rOPs were applied to detect heterogeneous vapour-liquid and liquid-crystal nucleation.

Simulation conditions
We truncated and shifted the LJ potential at r c = 2.5. van Meel et al. [49] report the coexistence data and vapourliquid surface tension for this model system in the relevant temperature range. As in [49], we fix the temperature at T = 0.45, in which case the vapour-liquid and vapourcrystal coexistence pressures are P coex vl = 4.28 × 10 −5 and P coex vx = 2.28 × 10 −5 , respectively, and the vapour-liquid surface tension is γ vl = 1.07.
We present results for two vapour pressures, P v = 1 × 10 −4 and P v = 3 × 10 −4 , with the corresponding bulk densities ρ L = 0.905 for the liquid and ρ S = 0.989 for the fcc solid [49]. For both the liquid and solid phases, Table 1 lists the associated vapour supersaturation and the difference in chemical potential to the vapour phase. At such low pressures, the free energy of the high-density phases does not change noticeably with pressure. Therefore, once the coexistence pressure is known, the difference in chemical potential can be computed directly from the vapour supersaturation.

Results
As was shown in [49,50], the bulk vapour-crystal nucleation for the present model system proceeds via Table 1. For the two vapour pressures used in this work, the vapour supersaturation S = P/P coex and the difference in chemical potential μ with respect to the vapour phase are presented for both the liquid (subscript l) and the fcc solid (subscript s) phases. The vapour-liquid coexistence pressure is P coex vl = 4.28 × 10 −5 and the vapour-solid coexistence pressure is P coex vs = 2.28 × 10 −5 . The difference in chemical potential between fcc solid and liquid is μ sl = μ s − μ l = −0.29. Subsequently, a crystal nucleus emerges within the liquid droplet. The two events happen on very different timescales and are, therefore, effectively decoupled except that the droplet needs to exceed a minimum size in order to sustain a stable crystal nucleus. As we show below, the pathway for heterogeneous nucleation follows essentially the same route as its homogeneous analogue, which means that the nucleation of the liquid and of the crystal are decoupled and can be treated independently.

Vapour-liquid nucleation
As a first step, we compute the vapour density as function of the distance perpendicular to a planar wall for various wallparticle interaction strengths . This computation is performed with a supersaturated vapour under the same conditions used in the subsequent nucleation studies. From this preliminary calculation, we can estimate the range of adsorption strengths for which heterogeneous liquid-vapour nucleation (and, a fortiori, liquid-crystal nucleation) is negligible on the timescale of our simulations. The results are presented in Figure 2 for a vapour pressure P v = 3 × 10 −4 (left) and P v = 1 × 10 −4 (right). The vapour density is rescaled by the bulk vapour density ρ 0 (P). As is to be expected for a dilute vapour, we find the vapour density to scale with the wall potential U w (R) according to the Boltzmann factor, ρ(R) = ρ 0 (P) exp {− β U w (R)} (data not shown). For stronger adsorption, where the vapour density is strongly increased, this expression requires corrections to capture the increasing inter-particle correlations. For wall-particle interaction strengths > 5 (P v = 3 × 10 −4 , left) and > 7 (P v = 1 × 10 −4 , right) spontaneous vapour-liquid nucleation was observed. Next, we compute the free-energy barriers for heterogeneous nucleation on a planar weakly attractive wall. It is important to note that the definition of a free-energy barrier for heterogeneous nucleation requires some care as, unlike the nucleation rate, it is not an experimental observable: the height of the free-energy barrier may depend on the choice of the reaction coordinate, in particular when different reaction coordinates are not linearly related (see e.g. [51]). However, if we choose the number of particles in the nucleus as reaction coordinate, then the free-energy barrier provides us with a measure of the probability to observe a nucleus of size n in a given system. This probability is made intensive by dividing the expected number of nuclei of size n by N, the number of particles in the system where nucleation takes place. If nucleation takes place in the bulk, N is simply the total number of particles in the system. However, if nucleation occurs on a wall, a result independent of system volume can be obtained by dividing the number of nuclei of size n that are formed on the wall by the total number of particles in contact Rep. ε = 1 ε = 3 ε = 5 ε = 7 Figure 2. (Colour online) Vapour density perpendicular to a weakly attractive planar wall for various wall-particle interaction strengths at a vapour pressure P v = 3 × 10 −4 (left) and P v = 1 × 10 −4 (right). For reference, a purely repulsive wall (Rep.) is shown. The repulsive potential has the same functional form as the attractive potential, but is truncated and shifted at the potential minimum to remove the attractive part.
with the wall. When comparing the rates of homogeneous and heterogeneous nucleation, one should then take into account how many particles in the system are in contact with the wall. In this way, one accounts for the fact that the rate of heterogeneous nucleation is proportional to the total 'active' surface area in the system (in other words, that it is proportional to the density of seeds). In what follows, we follow the above approach and define the barrier for heterogeneous nucleation as G het (n) = −k B T ln P het (n), where we define P het (n) as the ratio of M het (n), the number of clusters with size n in contact with the heterogeneity, divided by the (average) total number of particles in contact with the heterogeneity, N het . We identify a particle or cluster to be in contact with the heterogeneity using the restricted order parameter as described in Section 2.
The resulting free-energy barriers are presented in Figure 3 as function of droplet volume (to be more precise, we show the barrier as a function of the cube root of this volume. In what follows, we denote this quantity as the 'effective droplet radius'). The volume V is defined as the ratio of the number of particles n in the liquid nucleus, divided by the bulk density of the liquid ρ L = 0.905 (see Section 2 and [49]). Again, the vapour pressures are P v = 3 × 10 −4 and P v = 1 × 10 −4 for the left and right panels, respectively. When the attraction due to the wall is very weak, the computed barriers converge to the homogeneous case, as is to be expected. The barrier for homogeneous nucleation is plotted for the sake of comparison in Figure 3. With increasing attraction, the liquid-vapour nucleation barrier decreases until spontaneous nucleation occurs (not shown-as, in this case, it becomes impossible to measure a free-energy barrier). For the rest of this work, we choose our interaction strengths to be = 5 for P v = 3 × 10 −4 and = 7 for P v = 1 × 10 −4 . Outside the patches, the wall is purely repulsive.
Under these conditions, the size of the critical nucleus is still sufficiently small to be studied in US simulations for the system sizes that we used.
To obtain a rough estimate of the wall-liquid contact angle, we fit the expression from the classical nucleation theory (CNT) [52] with our simulation data (see Appendix 1 for details on the CNT expression), which yields = 0.33π (P v = 1 × 10 −4 , = 7) and = 0.43π (P v = 3 × 10 −4 , = 5). The CNT estimates for the heterogeneous nucleation barrier discussed below are based on these estimates of the contact angle. As our estimates of the contact angles are not very accurate, the comparison with our simulation results is mainly intended to gain a qualitative understanding of the observed trends. Several previous studies reveal that the CNT still hold for heterogenous nucleation on a planar wall if linetension effects are considered [53,54]. Thus, a quantitative comparison would require more precise knowledge of the contact angles and the solid-liquid-vapour line tensions.
Having established the limiting cases, i.e. homogeneous nucleation and heterogeneous nucleation on a planar wall, we focus on the effect of a planar circular weakly adsorbing patch on the nucleation rate. In particular, we are interested in the dependence of the nucleation rate on the radius of the patch. Figure 4 shows the nucleation barriers as function of the effective droplet radius V 1/3 . The theoretical CNT prediction is shown in the left panel and our simulation results in the right panel, both for a vapour pressure P V = 1 × 10 −4 . Valencia and Lipowsky [17] predicted (on the basis of CNT) that nucleation could be a two-step process with a metastable intermediate state under certain conditions. A similar conclusion was reached by Smorodin [16] on the basis of similar theoretical arguments and these intriguing observations stimulated recent theoretical explanation of contact line pinning effects for nanobubble stability [55]. Indeed, under the conditions studied in the present simulations, CNT predicts a two-step barrier for the vapour density P v = 1 × 10 −4 and a patch radius around R = 6.0, as depicted by the inset in left panel of Figure 4. However, our CNT results suggest that the metastable intermediate is likely to be found in a narrow range of patch radii that modulation of the barrier height is very small compared to the barrier height. In fact, the simulation results displayed in the same figure show no evidence for a two-step nucleation process. The dependence of the height of the nucleation barrier on the patch radius is presented in Figure 5. The nucleation barrier decreases monotonically from the homogeneous limit for a very small patch radius to the heterogeneous planar wall limit for large patch radii. The planar wall limit is reached at a patch radius R max ≈ R * homo sin , which corresponds to the contact circle radius of a critical nucleus on a planar wall. The CNT predictions are shown for reference. Qualitatively, the CNT predictions agree well with our simulations results.
Finally, we consider a pore geometry as shown in Figure 1, which is modelled as an open, weakly attractive hemispherical cavity embedded in an otherwise purely repulsive (i.e. hard) planar wall. Figure 6 presents the nucleation free-energy barriers as function of the effective droplet diameter for a vapour pressure P V = 1 × 10 −4 . Again, the left panel contains the CNT prediction, and the right panel the simulation results. As expected, the nucleation barriers approach the corresponding limiting cases for very large and very small pores. Interestingly, for intermediate pore sizes, CNT predicts a double barrier corresponding to a     two-step process, and this feature is indeed reproduced by our simulations. From a visual inspection of simulation snapshots of configurations at the local and global maximum of the free-energy barrier, we identify the first step to be a pore-filling event and the second step to be the nucleation out of the pore into the bulk. Figure 7 presents simulation snapshots along such a pathway. It shows that first a critical nucleus is required to nucleate a liquid inside the pore (left), the system then spends a certain time in the metastable filled pore state (centre), and finally, a second nucleation event is required to break out of the pore (right). This observation supports the predictions of Page and Sear for nucleation in a pore for an Ising model system [25]. The heights of the nucleation free-energy barrier are shown in Figure 8 as function of pore size, for both P v = 3 × 10 −4 (left) and P v = 1 × 10 −4 (right). The bottleneck for the overall nucleation rate is the crossing of the higher free-energy barrier shown in this plot. Again, good qualitative agreement with CNT is found. For nucleation in a pore, the dependence on the pore radius R is clearly nonmonotonic, featuring a global minimum of the free-energy barrier. This finding supports the suggestion by Chayen and co-workers [35] that, in a porous medium with a broad distribution of pore sizes, only a small number of pores (namely those that correspond to the minimum in the overall nucleation barrier) dominate the overall nucleation process.

Liquid-crystal nucleation
We now focus on the liquid-crystal nucleation. van Meel et al. and Chen et al. [49,50] showed that, in the case of homogeneous nucleation, the nucleation of a crystal in a previously nucleated droplet is relatively fast. In other words, the liquid-crystal nucleation event is not the rate-limiting  step. Moreover, van Meel et al. and Chen et al. [49,50] found that a droplet of a minimum size is required to host a crystal nucleus, and that this is approximately the size of a supercritical crystal nucleus plus a liquid mono-layer. Following the same procedure as in the case of vapourliquid nucleation, we first assess the limiting cases: homogeneous nucleation and heterogeneous nucleation on a planar wall. The crystallisation rate within a liquid droplet was reported in [49] and found to be independent of the droplet size and vapour pressure. For the sake of completeness, we here compute the crystallisation free-energy barrier for homogeneous nucleation in a bulk liquid in the constant number, pressure and temperature (NpT) ensemble. When we compare the present results with those reported in [49], we find good agreement.
For the crystal nucleation on a planar wall, we find that our system crystallises spontaneously on the timescale of our simulations for both attraction strengths, = 5 and = 7. Similar behaviour was observed in a study by Page and Sear on the pre-freezing transition of the same model system [56]. This high rate of spontaneous crystal nucleation implies that crystallisation may take place during US simulations that probe the free-energy barrier for vapour-liquid nucleation. To find out whether this is indeed the case, we analysed whether the configurations generated during the US simulations of the pathway for vapour-liquid nucleation contained any crystalline particles. In the left panel of Figure 9, we plot the number of crystalline particles n x , identified using bond-order parameters [6], versus the number of liquid particles n l , identified using Stillinger's overlapping spheres criterion [44]. Two observations can be made. First of all, all larger droplets (n l > 200) contained a crystalline cluster. If we consider that a crystal is always covered by a thin layer of liquid, as was found in [49], we can fit the results to an empirical formula assuming that the crystal takes the same shape as the liquid droplet, where c is a constant offset and the difference in density between the liquid and the crystal is ignored. The thickness of the liquid layer can then be extracted by dividing a by the geometrical prefactor linking the droplet volume to the radius of a spherical cap with liquid-wall contact angle , = a π 3 (1 − cos ) 2 (2 + cos ) −1/3 . Our fit yields a = 0.748 ± 0.1. Assuming a wall-liquid contact angle of = 0.33π , obtained from matching the simulation barrier height for vapour-liquid nucleation on a planar wall to the CNT prediction, we arrive at a thickness of ≈ 0.88, which agrees well with the assumption of a liquid mono-layer surrounding the crystallite. The difference a = n 1/3 l − n 1/3 x is also plotted in the right panel of Figure 9 and appears to converge to a constant value for large clusters.
Second, we observe that small droplets, n < 200, contain very little crystalline order. This means that for these small droplet sizes, crystallisation will not affect the computed vapour-liquid free-energy curve. As in the case of homogeneous nucleation, the liquid droplet has to exceed a critical size before crystals can form easily within it. In view of the high (spontaneous) crystal nucleation rate in large droplets, this might seem to be a rather surprising result, because the crystal critical nucleus is expected to be very small. However, due to the small wall-liquid contact angle, a droplet with n ≈ 200 is only a few layers thick. Crystal nucleation in such a pancake-shaped droplet is not determined by the total number of particles in the droplet, but by its thickness. Inside thin droplets, liquid-crystal nucleation is expected to be strongly suppressed, and that is indeed what we observe.
Next, we assess crystal nucleation in a finite-sized circular patch. Again, we observe spontaneous crystal nucleation for all but the smallest patch radii. Using US, we estimated the crystal nucleation barrier to be very small ( G * < 10k B T). Figure 10 shows the results from analysis of the crystallinity inside the liquid clusters. From the plot we conclude that, as in the case of nucleation on a planar wall, a minimum droplet size is required to support a stable crystal. As before, we find that the resulting crystal fills the droplet except for a liquid mono-layer. The fact that spontaneous crystal nucleation in larger droplets is observed even for the smallest patch radii (i.e. R = 3) suggests that the structure of the substrate strongly influences crystal nucleation.
For crystallisation in a pore, the curvature can cause stress and frustration for a regular arrangement of particles, and thus might inhibit crystal nucleation. Therefore, we computed the free-energy barrier for crystal nucleation for those systems within a pore filled with a metastable droplet (data not shown). For small pores, no stable crystal grows in these droplets. For the largest pore that can contain metastable liquid droplets (R = 6), we observe a low crystal nucleation barrier of G * ≈ 7.5k B T. This barrier height is slightly less than the nucleation barrier for growing out of the pore ( G * ≈ 10k B T). To identify whether crystallisation precedes pore break-out, we need to compare transition rates rather than free-energy barriers. Using a diffusive kinetic prefactor, which reflects the diffusive nature of our Monte Carlo simulations, we can express the rate according to [49]: Here D is the self-diffusion coefficient, ρ the density of the parent phase and N refers to the number of monomers that can initiate the transition. Note that in the case of crystallisation, N equals the number of particles in the liquid phase, whereas in the case of pore break-out, the droplet counts as a single 'monomer' that can initiate the transition; hence N = 1. The diffusion constant can be approximated by the simulation step size squared per Monte Carlo cycle, D ≈ x 2 /τ . This leads to a pore break-out rate for the liquid of J l ≈ 4 × 10 −11 τ −1 and a crystallisation rate of J c ≈ 4 × 10 −5 τ −1 . Hence in these pores, crystallisation is expected to precede pore break-out. This prediction is confirmed by direct simulations. The resulting stable crystal contains approximately n = 140 particles.
For pore radii R > 6, a fully filled pore is not even metastable and the liquid continuously grows into the bulk. In this situation, the pore's attraction and curvature can still influence the rate for crystal nucleation, but this step is never rate-limiting.

Conclusion
The focus of the present paper is on heterogeneous nucleation of liquids and crystals. We studied nucleation on a planar wall, on a planar circular patch and in a hemispherical pore. For all surface geometries, we find the nucleation pathway to follow the homogeneous two-step route, with the vapour-liquid nucleation as rate-limiting step. Comparing the simulation results to the corresponding CNT predictions, we find good qualitative agreement. However, a truly quantitative comparison was not possible due to the fact that the wall-liquid and wall-crystal surface tensions are not known with sufficient accuracy, and the line tensions of the nuclei on the surface are unknown. Yet line tensions are known to have a pronounced effect on the nucleation barrier [8,53,57,58].
For a planar circular patch, the free-energy barrier for vapour-liquid nucleation shows a monotonic decrease with increasing patch radius, interpolating between the barrier for homogeneous nucleation and the one for nucleation on a planar wall. For the parameter range studied in our simulations, we did not observe two-step vapour-liquid nucleation. Moreover, CNT calculations using as input, the parameters that characterise the systems that we studied in the simulations also failed to exhibit two-step nucleation. These findings indicate that the two-step nucleation scenario on a planar circular patch that was predicted in [16,17] occurs for conditions that are outside the domain that we could study in our simulations. Furthermore, our simulations show that liquid-crystal nucleation happens spontaneously for large patches and at a planar wall, provided that the liquid droplet exceeds a minimum size of n ≈ 200 particles.
For a hemispherical pore we find a richer behaviour. Very small and very large pores exhibit a vapour-liquid nucleation free-energy barrier that approaches the ones for homogeneous nucleation and heterogeneous nucleation on a planar wall, respectively. For intermediately sized pores, the liquid nucleation itself is a two-step process. In the first nucleation event, a droplet nucleates inside the pore and grows to fill the entire pore. From this metastable state, a second nucleation event is required to grow out of the pore. The overall nucleation rate shows a non-monotonic dependence on the pore radius with an optimum pore size for which the barrier is the lowest. This observation corroborates a previous study on nucleation in and out of pores for an Ising model system [25]. Depending on the strength of attraction, the minimum barrier can be even lower than for nucleation on a planar wall. The liquid-crystal nucleation typically follows once the droplet is large enough. For a small pore, the droplet extends into the bulk and the crystal nucleation is homogeneous. For large pores, the crystal forms in contact with the pore surface at a very high rate -in practice, there is no nucleation barrier. Only for the largest pore within which a metastable droplet could be maintained, we did find the filled pore to be large enough to host a stable crystallite. In that case, the liquid crystallises before it grows out of the pore.
In conclusion, our simulations show that a weakly attractive circular patch or hemispherical pore on an otherwise non-adsorbing substrate can reduce significantly the free-energy barrier for nucleation. Both geometries are most efficient if their radius is close to or slightly larger than the critical radius for homogeneous nucleation, at which point they strongly accelerate the formation of exactly one nucleus. A widely spaced array of adsorbing patches or pores with a broad distribution of radii crafted on an otherwise non-adsorbing surface should, therefore, make an efficient nucleation agent for the formation of high-quality crystals. We point out that the present work is limited to smooth surfaces. Further research is required on (crystal) nucleation at surfaces that are rough on the scale of a particle diameter. In such cases, roughness is expected to have a pronounced effect on the free-energy barrier for crystal nucleation. Figure A1. Left: liquid droplet on a planar surface, surrounded by vapour. The droplet's shape corresponds to a spherical cap where the geometrical parameters are the droplet radius R, the contact angle , the height of the droplet h and the radius a of the contact circle between droplet and surface. Right: liquid droplet in a spherical pore, surrounded by vapour. The droplet is composed of two spherical caps with the corresponding parameters R i , a i , h i and i , with i = 1 for the lower cap and i = 2 for the upper cap. The contact angle is denoted (without subscript). For the relation between these parameters, see the main text.
interfaces by σ wg , σ wl and σ lg , respectively, and the contact angle by . The CNT expression for the free energy of this droplet is given by [1] G = σ lg A lg − | μ|ρV l + σ wl − σ wg A wl , where A lg and A wl refer to the vapour-liquid and wall-liquid surface areas, respectively, V l to the droplet volume, ρ to the bulk liquid density and μ to the difference in chemical potential between the bulk vapour and bulk liquid. Using the Young-Dupré relation for the contact angle, σ wg = σ wl + σ lg cos , we obtain G = σ lg A lg − | μ|ρV l − σ lg A wl cos (A1) The volume and the surface area of this droplet are given by the expression for a spherical cap: where V s (R) = 4π /3R 3 and A s (R) = 4π R 2 is the volume and surface area, respectively, of a full sphere with radius R, and V ( ) = 0.25(1 − cos ) 2 (2 + cos ) and A ( ) = 0.5(1 − cos ) are contact-angle-dependent volume and surface modifiers, respectively. Inserting the volume and area of our spherical cap in the CNT expression yields the barrier height and critical radius for nucleation on a planar attractive wall [1]: with G * homo = 16π/3σ 3 lg ρ −2 | μ| −2 the critical CNT barrier height for homogeneous nucleation. The critical nucleus radius is the same for both homogeneous nucleation and nucleation on a planar wall, R * = 2γ ρ −1 | μ| −1 , and the radius of the wall-liquid contact circle is R = R * sin .

Circular patch
For an attractive circular patch with radius R P , the free energy has to be slightly modified. We can identify three distinct regions [16,17,59]: in regime I (R < R P /sin P ), the droplet is small enough to fit entirely on the circular patch with the natural patch-liquid contact angle P . Regime II applies if the droplet covers the entire patch with a contact angle P < < W , where W is the wallliquid contact angle. In the final regime III (R > R P /sin W ), the droplet extends beyond the circular patch, and the contact angle is that with the wall. The CNT expression for the free energy difference is hence G I (R) = σ lg A S (R) A ( P ) −σ lg πR 2 sin 2 P cos P −| μ|ρV S (R) V ( P ) G II ( ) = σ lg A S (R P / sin ) A ( ) −σ lg πR 2 P cos P −| μ|ρV S (R P / sin P ) V ( ) G III (R) = σ lg A S (R) A ( W ) −σ lg πR 2 P cos P −σ lg π (R sin W − R P ) 2 cos W −| μ|ρV S (R) V ( W )

Hemispherical pore
Next, consider a droplet in a hemispherical pore, as sketched in the right panel of Figure A1. The droplet volume consists of two spherical caps, each with its own radius R i and contact angle i . From the figure, it is obvious that P = 1 + 2 , where P is the liquid-pore contact angle. The lower cap has the curvature of the spherical pore, so R 1 = R P . Both caps share the radius of the circle where they meet, so a 1 = a 2 ≡ a. The radius R i of a cap is related to a by a = R i sin i , so the radii are related by R 2 = R P sin 1 /sin 2 = R P sin 1 /sin P − 1 . Analogous to the droplet on a circular patch, we identify three regimes: in regime I, the droplet does not fill the pore, i.e. 1 < π/2. Note that this is independent of the actual contact angle P . In the second regime, the droplet fills the pore, 1 = π /2, and grows until the upper cap angle 2 reaches the contact angle with the surrounding wall, P − π /2 ≤ 2 ≤ W . Here, we point out that for intermediately to strongly wetting pores, P < π/2, the upper cap angle 2 = P − π /2 starts out negative. In regime III, the upper cap grows beyond the pore onto the surrounding wall, taking a mushroom-like shape, R > R P . Computing the free energy for the three different regimes then yields G I ( 1 ) = σ lg A S (R P ) A ( P − 1 ) sin 2 1 sin 2 P − 1 −σ lg A S (R P ) A ( 1 ) cos P −| μ|ρV S (R P ) V ( P − 1 ) sin 3 1 sin 3 P − 1 −| μ|ρV S (R P ) V ( 1 ) G II ( 2 ) = σ lg A S (R P ) A ( 2 ) sin −2 2 −σ lg A S (R P ) A ( 1 ) cos P −| μ|ρV S (R P ) V ( 2 ) sin −3 2 −| μ|ρV S (R P ) V ( 1 )