# Computer simulation studies of wetting on heterogeneous surfaces

## I Introduction

The wetting of solid surfaces by fluids is a problem of great practical importance that has been extensively studied over the years . Most often, the experimental work has involved measurements of the contact angle $`\theta `$ made by a liquid on the solid surface of interest. Macroscopically, this is the angle made by the surface of a droplet of liquid in contact with a solid surface and in equilibrium with its vapor. Young’s equation gives the relation between $`\mathrm{cos}\theta `$ and the interfacial tensions $`\gamma _{LV}`$, $`\gamma _{SL}`$ and $`\gamma _{SV}`$ where L, S and V denote liquid, solid and vapor respectively:

$$\mathrm{cos}\theta =\frac{\gamma _{SV}\gamma _{SL}}{\gamma _{LV}}$$
(1)   
(Microscopic expressions relate these surface tensions to spreading pressure, which is $`p_{}`$, the component of the pressure tensor parallel to the interface.) The point here is that when the ratio on the right hand side of eq. (1) is larger than one, a droplet will spread into a thin layer (and the equation actually becomes an inequality, since $`\mathrm{cos}\theta `$ cannot be greater than 1 or less than -1) and when the ratio is less than -1, the droplet will not spread. These two regimes are described as wetting and non-wetting respectively and the regime where $`\mathrm{cos}\theta `$ is finite is called partial wetting. Note that the molecular equation that gives interfacial tension as an integral over the z-dependent local $`p_{}`$ is well established only for those systems where the fluid-solid interaction is independent of the variables $`x,y`$ that give the position of the fluid molecule over the solid. This limitation is a major obstacle to the use of this formalism in microscopic calculations of the wetting properties of fluids on rough or heterogeneous surfaces.

However, many molecular simulations that conform to the flat surface limitation have been reported . An example  is shown in Figure 1, where values of the contact angle are shown for a Lennard-Jones fluid (i.e., one with potentials that depend upon the inverse 12,6 powers of separation distance) interacting with a surface with a potential that varies as $`z^9`$ for its repulsive part and $`z^3`$ for its attraction. Results are shown for two values of temperature $`kT/ϵ`$, where $`ϵ`$ is the interaction energy well depth for the fluid-fluid interaction. The quantity $`ϵ_W`$ is one of great importance in the present discussion since it gives the ratio of the fluid-solid well depth to that for the fluid-fluid interaction. It turns out that the values of this parameter are generally the most important factor in determining wetting behavior. Note that the two curves in Figure 1 indicate the ranges of $`ϵ_W`$ that give partial wetting, with complete wetting and drying occurring for $`ϵ_W`$ larger and smaller, respectively, than the interval where $`\mathrm{cos}\theta `$ is between 1 and -1.

An alternative to the contact angle computation is the more generally applicable formulation of this problem based on simulations (and measurements) of adsorption isotherms for a fluid on a weakly interacting solid surface. In this case, the limits of wetting and non-wetting are easily stated: wetting is the normal case for a strongly interacting surface where the adsorbed amount increases smoothly to infinity as the pressure approaches the bulk vapor pressure of the fluid. The non-wetting case is also easily described, since the adsorption isotherms in this case never rise above a rather small level (usually less than a monolayer) until they meet the bulk vapor pressure line where the container will fill with bulk fluid. An example  is shown in Figure 2, which gives simulated isotherms for neon adsorbed at several temperatures on a weakly interacting surface. All of these show a vertical rise at the bulk vapor pressure starting from a small coverage typical of non-wetting behavior. The estimated monolayer capacity for neon on this surface is 9 atoms/nm<sup>2</sup> . (The vertical rises terminate at finite values because of the finite capacity of the box used in the simulations.) The temperature range for these isotherms is from somewhat higher than the bulk boiling point ($`27.1`$ K) to close to critical ($`44.4`$ K).

The analogue to partial wetting in systems with finite contact angles is not as clear as for the wetting and non-wetting limits. In adsorption, one has the phenomenon known as prewetting which appears to be the adsorption version of partial wetting. In prewetting, the adsorption isotherm makes a finite jump in coverage at some value of the pressure less than the bulk vapor pressure $`p_o`$. At coverages above this transition, the isotherm is close to that for a wetting system, with coverages approaching infinity as $`pp_o`$. Below the transition, the isotherm has the general appearance of one for a non-wetting fluid. Examples  are shown in Figure 3, which also shows a single non-wetting isotherm. All are for the temperature of $`28`$ K corresponding to a bulk vapor pressure of $`1.23`$ atm.; the monolayer capacity is $``$ $`1.4`$ particles per unit cell.

Although the quantitative connection between prewetting isotherms and partial wetting contact angles has not been established yet, it is evident that this need not prevent one from using adsorption isotherms to characterize the wetting behavior of fluids on weakly interacting surfaces. After giving a brief description of the computer simulation algorithms used to evaluate the thermodynamic properties of adsorbed films, the remainder of this paper will be devoted to a presentation of results for simple Lennard-Jones atoms (rare gases) adsorbed on a number of model weakly interacting solid surfaces, some of which are heterogeneous.

## II Computer Simulation Algorithms

Computer simulation of the macroscopic properties of molecular fluids is now a very widely used technique to solve the problem of evaluating the microscopic equations for these properties. These simulations fall into two general categories: first, molecular dynamics, where the equations of motion of the molecules are solved numerically to generate their configurations as a function of time. This data can then be used to evaluate both thermodynamic and dynamical properties of the system. Second, one can utilize Monte Carlo algorithms which also generate configurations using a probabilistic approach. The Monte Carlo algorithms have proved to be very effective in problems relating to physical adsorption. Therefore the discussion here will be limited to such computations. We need only briefly recapitulate the general principles, since these calculations have been extensively described in monographs  and review articles.

In Monte Carlo simulations, atoms are subjected to random moves, which are accepted or rejected according to rules that will, in the limit of a large number of trials, produce configurations that correspond to those in the desired statistical mechanical ensemble. We discuss briefly here the requirements for generating a canonical ensemble, and then the more useful Grand Canonical ensemble. A canonical ensemble consists of a collection of $`N`$ particles in a box of volume $`V`$ at temperature $`T`$. Interaction potentials between particles are chosen so that one can evaluate the total potential energy of the entire system of particles. Usually, pair-wise interactions are specified and for atoms, the energy depends only upon their position coordinates in the box. The system is extended by using periodic boundary conditions and the minimum image convention, which means that the computer box of several hundred or thousand particles is surrounded by images of itself and the box walls are made transparent so that atoms can pass out (and images pass in). Each atom interacts with its nearby real neighbors as well as with nearby images. The crucial question is: when does one accept (or reject) a move than produces a change $`\delta U`$ in the potential energy of the system? The answer is in two parts: If $`\delta U`$ is negative, the system is more stable and the move should be accepted; if it is positive, then it may or may not be accepted depending upon the value of $`\mathrm{exp}(\delta U/kT)`$. If this is larger than a randomly chosen number between 1 and 0, the move will be accepted. Otherwise, it is rejected and another random move is tried. After a large number of such trials, the configuration of the atoms in the system will conform to that of the atoms in a canonical ensemble and one can use such configurations to evaluate averages of quantities such as $`<U>`$, the mean potential energy. Note that the mean kinetic energy is $`3NkT/2`$ in such a system. The average pressure can also be evaluated using the virial theorem; more to the point, the elements of the pressure tensor parallel and perpendicular to the surface in an inhomogeneous fluid such as an adsorbed layer can also be evaluated.

This brief discussion of a canonical ensemble Monte Carlo has been given primarily to help establish the ground rules for the Grand Canonical ensemble Monte Carlo (GCMC) of a physically adsorbed gas in which the quantities to be evaluated are $`<U>`$ and $`<N>`$, the average potential energy of the adsorbed fluid and the average number of particles in a system with fixed volume $`V`$, temperature $`T`$ and chemical potential $`\mu `$. For an adsorbed layer in equilibrium with gas at pressure $`p`$, the chemical potentials of the two phases must be equal and, assuming that the gas phase is ideal (or that corrections to $`\mu _{gas}`$ for non-ideality of the gas phase can be made), one has

$$\mu _{gas}=\mu _{ads}=kT\mathrm{log}p+q$$
(2)   
where q is a constant calculable from the standard expression for an atomic gas. Thus, the GCMC simulation treats $`<N>`$ and $`<U>`$ as quantities to be evaluated from the simulation at fixed $`V`$, $`T`$, $`\mu _{ads}`$ , with $`\mu _{ads}`$ related to the observable pressure by eq. 2. The moves in this Monte Carlo procedure are of three kinds: translation, as in the canonical case; creation of a new atom at a randomly chosen position in the computer box; and deletion of one of the pre-existing atoms. The decision concerning the acceptability or not of a given move now depends upon the value of $`\delta C`$ or $`\delta D`$ rather than $`\delta U`$. $`C`$ and $`D`$ refer to creation or destruction, respectively and are equal to:

$$\delta C=\delta U+\mathrm{log}(zV/(N+1))$$
(3)   
$$\delta D=\delta U+\mathrm{log}(N/zV)$$
(4)   
where N is the instantaneous number of particles in the system and z is the activity defined by $`\mathrm{exp}(\mu /kT)/\mathrm{\Lambda }^3`$ , where $`\mathrm{\Lambda }=h/\sqrt{2\pi mkT}`$. (For an ideal gas, $`z=p/kT`$.) By combining moves of these three types, one can generate configurations in the Grand ensemble after a considerably larger number of moves than for the analogous canonical ensemble. At this point, the isotherm is essentially a plot of $`<N>`$ versus $`kT\mathrm{log}p=\mu _{ads}+`$ constant.

This formalism is very widely used in simulations of adsorption, including those for weakly interacting solid where wetting is an important issue. The wetting of polar molecules on non-polar adsorbing surfaces has been studied for many years. The simpler case of rare gases on real weakly adsorbing surfaces did not receive much attention until it was pointed out that the alkali metals should have gas-solid interactions that are sufficiently weak to give interesting wetting behavior. The well depth ratios for several of these systems are listed in Table I. Both experiment and theory indicate that many interesting features are exhibited by these systems.

The remainder of this paper will be devoted to a description of some of the results obtained for the wetting behavior of neon adsorbed on these weakly interacting solid surfaces. Several simulations have been reported for a solid called CO<sub>2</sub> \- a ratio is listed for this material as well, although one should not take its description as solid CO<sub>2</sub> seriously. From Table I, it is evident that the alkali metals Cs and Rb are expected to be the weakest possible adsorbents for rare gases. A number of wetting studies have been carried out for Neon on these materials , primarily because one hopes that neon is sufficiently heavy to obey classical statistical mechanics, which is the basis for the GCMC algorithms used.

## III Simulation Results

Simulations were first performed for the energetically flat approximations to the gas-solid interactions because of their simplicity. Figures 2 and 3 already show some simulations of neon on model alkali metal surfaces. Note that the simulations were not restricted to interaction well depths predicted by theory. Trends in the adsorption behavior could most easily be observed by taking various interaction strengths that do not necessarily correspond to the theoretical predictions for a given metal.

Simulations for Ne on a surface with $`ϵ_W`$ = 0.5 are shown in Fig. 2 and indicate that this system is similar to the Cs case of Table I in that they also are non-wetting over a range of temperature up to a couple of degrees below critical. (Simulations closer to critical than these were not feasible for the usual reason: correlation lengths were beginning to be larger than the computer box and were therefore introducing artifacts in the isotherms.) We note that one can readily extract local densities in the adsorption box from the configurations of the atoms in it. In this work, a box of lateral dimensions = $`2.782.78`$ nm<sup>2</sup> and height 7.5 nm was used. In the lateral directions, periodic boundary conditions were applied; in the ”vertical” or $`Z`$ direction, the adsorbing surface was a wall at $`Z`$=0 and a hard wall was imposed on the opposite side of this box. The local densities are shown in Figure 4 for a non-wetting case, in Figure 5 for a wetting example and in Figure 6 for a prewetting system. For a non-wetting system, the figure shows that the densities of atoms at both the lower and upper walls is very small, as one might expect. (This behavior has been observed previously by Adams and Henderson ). Thus, no complete monolayer is formed, giving very small adsorption over the entire pressure range up to $`p_o`$. If the gas-solid well depth is increased, one observes prewetting, as in Figure 6. The densities for this prewetting case are at least qualitatively similar to those for the wetting system shown in Fig. 5. Note that this prewetting transition involves a jump in coverage from about $`2.2`$ atoms per nm<sup>2</sup> (0.24 layer) to 30.3 atoms per nm<sup>2</sup> ($`3.3`$. layers) for $`ϵ_W`$=2.8 at T=28 K.The prewetting regime for a solid with well depth of $`95`$ K ranges from $``$ $`22`$ K to $``$ $`29`$ K.These results appear to depend somewhat upon the width of the potential well in addition to its depth .

As for the explicitly heterogeneous surface, we start from the assumption that any variation in gas-solid energy across the surface can be described as heterogeneity. Thus the periodic energy variations due to the effect of atomic structure when the pair-wise additive approximation for the adsorption potential replaces the integrated perfectly flat surface falls into this category and thus were briefly considered . Figure 7 shows that the isotherm obtained for a periodic surface with average interaction energy equal to that

for the flat surface is noticeably different from that for the flat surface. What has happened is that the Boltzmann-weighted average gas-solid interaction for the periodic case is somewhat more negative than that for the flat surface, giving a Henry’s Law low coverage isotherm that is slightly steeper for the periodic case than for the flat. In particular, the prewetting jumps in Fig. 7 start from coverages of $`0.35`$ monolayers (periodic) and $`0.24`$ monolayers (flat) and end essentially at $`3.0`$ monolayers.

## IV Conclusions

The two approaches to wetting based on Young’s equation for the surface tensions and adsorption isotherms should be mutually consistent. It is likely that the boundaries for the surface tension wetting$``$partial wetting and partial wetting$``$non-wetting are indeed consistent with the isotherm results, if one assumed that prewetting is in fact an alternative notation for partial wetting. (Note also that non-wetting here means only slight adsorption before the jump to bulk phase condensation occurs, which has elsewhere been described as partial wetting, thus confusing the issue.) The points of resemblance between partial wetting from Young’s equation and prewetting from the isotherms have not been elucidated. Indeed, no detailed comparison of the wetting and non-wetting boundaries as a function of $`T^{}`$ and $`ϵ_W`$ for the two methods has yet been presented.

Not many simulations of isotherms showing prewetting have been presented until fairly recently. This is at least in part because the prewetting jumps appear to occur at pressure very close to the bulk condensation value. This not need be the case, since it has been shown that prewetting isotherms on surfaces with some roughness, either due to periodic variation in energy across the surface (see Figure 7) or to heterogeneity introducing by adding or deleting atoms from the surface to form defects will cause the prewetting jump to shift down to a pressure where the transition is more easily observed. However, for sufficient energy defects in the surface, the transition is no longer sharp; i.e., not first order. Still, the observation of truly first order transitions in computer simulations has always been difficult The situation here where one has no hint of where the first order behavior does occur other than from the simulated isotherms (or from the shape of $`\mathrm{cos}\theta `$ curves such as those shown in Fig. 1) is evidently one where the difficulty of determining first order transitions is a major one.