Impact of Surface Roughness on Crystal Nucleation

: 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 ﬂuctuations (rough-ness) 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.


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 threedimensional 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.

Simulation Model
In our model, the simulation cell is partitioned into a lattice of N = N x × N y × 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 × 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 −∆E/k B T ) where ∆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(t) 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, ∆t t = ν 0 /N t where the "attempt frequency" ν 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 ν 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 ν 0 so as to give the observed single-molecule diffusion constant. However, for simulations, one takes a to be the length scale and ν 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.

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(r) = 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) ≡ 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(h)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 σ, 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(r 1 , r 2 ) = h(r 1 )h(r 2 ) . 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 r 1 and r 2 , e.g., G(r 1 , r 2 ) = G(|r 1 − r 2 |) = G(r 12 ). In this case, a measure of the typical length scale of such correlations, i.e., of the lateral correlation length, is given by ∞ 0 G(r)r dr/ ∞ 0 G(r) dr. In our work, we have used the standard exponential model, where it can be confirmed that the parameter ξ 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 and in the high roughness limit, one gets

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, λ, 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 λ ≤ λ S where λ 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 λ X ≥ λ C for some (crystal) threshold λ 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 λ < λ S and end up in λ > λ C and is estimated as the product of Φ 0 , the initial flux of trajectories that cross λ S and of the probability that a system that has λ = λ S evolves into one with λ = λ 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 λ < λ S to λ = λ 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 λ 1 > λ S is defined and the probability of the transition λ S → λ 1 is determined by randomly choosing a configuration from the database of configurations with λ = λ S and evolving it until it either reaches λ 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 λ 1 gives the probability of the transition. This is then repeated for a series of such milestones, λ S < λ 1 < λ 2 < · · · < λ N = λ 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 λ S = 5, λ C = 100 and milestones at each integer in between. From each milestone λ k , new trajectories are launched until there are 400 that reach the next milestone.

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 ρ = 0.0016 in lattice units. The bulk density was chosen sufficiently low that the nucleation was not instantaneousmaking 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 σ 2 and correlation length ξ, typical examples of which are shown in Figure 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 ξ = 5 (upper graphic), and ξ = 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 σ 2 = 0.5 to σ 2 = 4.5. Moreover, at fixed variance the nucleation rate increases by a factor of a thousand as the correlation length decreases from ξ = 20 to ξ = 5. The critical size decreases almost monotonically for the smaller correlation length but at larger correlation length it is unaffected until the roughness exceeds σ 2 = 4. Our goal in the following is to understand the mechanisms behind these variations.

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] where γ sv , γ sc , and γ cv , are the interfacial tension for substrate-vapor, substrate-crystal, and crystal-vapor interfaces. By combining (5) with Young's equation γ sv − γ sc = γ cv cos θ 0 one obtains the Young-Dupre equation where θ 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 (r) = γ cv (1 − cos θ(r) ), where θ(r) is the contact angle modified by the local slope at the surface, θ(r) = θ 0 + arctan(|∇h|). 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 This shows that for fixed correlation length, adhesion energy increases with increasing variance, σ. 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 ρ s as a function of σ and ξ. The surface density ρ s is defined as the number of molecules whose distance from the plane z = 0 is lower so that h max = max({h(r) : r = (x, y)}) is normalized by the surface area of the plane, N x × N y . Figure 5 shows that ρ s is indeed proportional to σ, 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.

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 σ 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 λ = 100, as a function of σ 2 and ξ.
The probability distribution of A sc , obtained from a thousand substrate-cluster configurations, is displayed Figure 8. One observes that for σ 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 σ 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 σ 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(A sc ) 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 λ 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 | |/2a 2 whereas the crystal-substrate surface tension is γ 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 (σ 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 | ≤ 0.3.

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 substratedominated 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. Acknowledgments: The authors are grateful for discussions, insights and guidance from Pierre Gaspard. Computational resources have been provided by the Shared ICT Services Centre, Université Libre de Bruxelles.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Model of Rough Surface
Appendix A.1. Generating the Surface A homogeneous rough surface can be considered to be a two-dimensional random process with heights specified by a height auto-correlation function, G(r), or, equivalently, a spectral density, S(k) defined as This is the starting point of the Fourier transform method [18] used for the generation of random rough surfaces. Consider first a discrete one dimensional surface h n , n = 0, 1, ..., N − 1 whose Fourier transform givesh k , k = 0, 1, ..., N − 1. Specifically, define the surface as where the phases φ k are uniformly distributed between 0 and 2π and with the requirement that φ N−k = −φ k . Then, it follows from Parceval's theorem that ∑ N−1 n=0 h n h n+r = N ∑ N−1 k=0S k = S r so that one can reverse this and begin by specifying S r = G(r), then determining its Fourier transform,S k and then constructing h n which has the desired correlation function. Furthermore, since the heights are identically distributed, independent variables the central limit theorem assures that, in the limit of large N, the surface-averaged height h ≡ ∑ N−1 n=0 will be a Gaussian process with mean zero and variance S 0 = G(0). Similarly, for the case of a two-dimensions, a random rough surface h nm , with prescribed spectral coefficients S kl can be constructed fromh kl = √ S kl exp(iφ kl ), where as before φ kl are random phases. Then defines a Gaussian random surface with surface height distribution This means that for any h nm , the probability that the height is between h and h + dh is p(h) dh. At any position on the surface h(r) = 0, where · means an average over an ensemble of random surfaces, and therefore the variance σ 2 = h 2 characterizes the roughness in the direction normal to the surface. In the following we will assume the surfaces are ergodic random fields, meaning that all statistics can be obtained from a single large random surface.
For a homogeneous, isotropic surface, as assumed in this work, we expect to start with a model for S surf (k) = S( k 2 x + k 2 y ). However, in order to generate the surface using the expression above, we need the spectral densities S(k x ). These are related by an Abel transform [25,26] and the reconstruction of a circularly symmetric two-dimensional function from its projection is known as Abel inversion of the projection which is The morphology of an isotropic random can often be characterized by a single parameter, the lateral correlation length ξ = ∞ 0 G(r)r dr/ ∞ 0 G(r) dr. Using this, our Gaussian surfaces (see Equations (1) and (2)), are generated with a surface spectral density given by (A7) Since numerical random surface are of finite extend the discrete spectral density S kl in (A3) is a function of the system size N. However, for large system size, S k,l S surf (|k|). This is shown Figure A1 where S k,l=0 is displayed for N = 4096 together with S surf (k). where |∇h| is the local surface slope. For weak roughness, |∇h(x, y)| 1, 1 + |∇h(x, y)| 2 1 + 1 2 |∇h(x, y)| 2 − 1 8 |∇h(x, y)| 4 + · · · which yields In our model, the solid substrate consists of a stack of cubes whose upper surface is a collection plaquettes. The area of this discrete surface is given by where a 2 is the area of a plaquette and the height h is dimensionless. We shall qualify this expression as the strong roughness limit. Most of the roughness effects can be related to the average of the surface area (A9), which is an expansion in (2n, 2m)-order moments (∇ x h) 2n (∇ y h) 2m . If the surface height obeys the Gaussian distribution (A4), then the joint distribution of two heights h 1 and h 2 is where matrix M is the inverse of the matrix h i h j [27]. From p(h 1 , h 2 ) the following joint distribution of h and ∇h can be derived, [27,28] p(h, ∇h) = p(h) and, more specifically, where ρ = (∇h) 2 1 2 is the rms of the local slope. Using (A14) one obtains the following ensemble average, which, when inserted into (A11), gives Therefore the rms local slope alone determines the surface area averaged over roughness configurations. In order to obtain ρ 2 one introduces ∆h(r) = h(r 0 + r) − h(r 0 ) which, in terms of the the Fourier transform of the surface heighth(k), is ∆h(r) = h (k) e ik·r 0 (e ik·r − 1)dk. Then, [∆h(r)] 2 = A 0 e i(k+k )·r 0 h (k)h(k ) (e ik·r − 1)(e ik ·r − 1)dk dk , where A 0 comes from a sum over the position r 0 . Since, using translational invariance, giving [∆h(r)] 2 = S surf (k)(2 − e ik.r − e −ik.r )dk (A20) = ∞ 0 dk 2π 0 dθ k S surf (k)(2 − e ik r cos θ − e −ik r cos θ ) where we used the relation 2π ∞ 0 S sur f kdk = σ 2 and J 0 (x) = J 0 (−x). Finally, introducing the lattice spacing a in the definition of the space derivative, ρ 2 = lim r→a 2 a 2 [∆h(a)] 2 gives, using expression (2) for G(r), For a ξ this expression yields We generated 10 3 random rough surfaces of area A 0 = (256) 2 to obtained "exact" average quantities for comparison with analytical expressions. The results are displayed as functions of σ for two values of the correlation length (ξ = 5 and ξ = 20). Figure A2a displays the results for ρ 2 as given by (A22). From relation (A22) one obtains the following expansion for the average Wenzel factor as a function of σ, with b = 4/a. In Figure A2b the approximated excess surface area r exc = r w − 1 is displayed as a function of σ. Only the first three terms of the expansion (A23) are included in the analytical expression. In order to compare using the same "continuous" treatment of the surface, the figure compares this result to simulations of the surface from which the area was computed via linear interpolation from the center of one cell to the centers of the neighboring cells.
In the high roughness limit and treating the surface as being made up of square plaquettes as it actually is in the simulations, one takes the average of (A10) over the distribution (A14) which gives Inserting (A21) with a = 1 in this expression yields the high roughness average Wenzel factor (A25) Figure A2c shows a comparison of this analytic expression to the area as determined from simulating the surface and computing the area based on the cubic structure of the actual simulations. Figure A2. (a) The rms local slope, ρ, as predicted by Equation (A22) as functions of variance, σ, for the cases of ξ = 5 (black curves and symbols) and ξ = 20 (red). (b) A comparison between the analytic approximations of Equation (A23) (first three terms, solid lines) and the "exact" values of the excess relative surface area, r exc ≡ r w − 1 as functions of σ with the same labeling. The "exact" values were obtained by averaging over a thousand realizations of the surface and the area was computed in a "continuous" approximation as explained in the main text. (c) The same thing but treating the surface as being made up of discrete plaquettes for both the analytic approximation, Equation (A25), and in computing the area of the simulated surface.