J. Phys. Chem. C 2008, 112, 12905–12913
12905
Influence of Substrate Strength on Wetting Behavior Michael S. Sellers and Jeffrey R. Errington* Department of Chemical and Biological Engineering, UniVersity at Buffalo, The State UniVersity of New York, Buffalo, New York ReceiVed: April 21, 2008; ReVised Manuscript ReceiVed: June 7, 2008
We examine the evolution of prewetting phase behavior with substrate strength for a model Lennard-Jones system, which consists of monatomic particles interacting with a single structureless surface. Grand canonical transition matrix Monte Carlo simulation is used to characterize surface phase behavior for several wall strengths. Our results indicate that both wetting and prewetting critical temperatures increase with decreasing wall strength, with the former increasing at a faster rate than the latter. The two curves appear to meet at the critical temperature of the bulk fluid at nonzero substrate strength. Overall, prewetting transitions are found within a triangular-shaped region within the temperature-wall strength plane, bound by wetting and prewetting critical lines and a series of layering/crystallization transitions at low temperature. The interfacial tension is calculated along each saturation line, from which we deduce the latent heat of the prewetting transition. Several corresponding states ideas are also discussed. We examine the ability of the phenomenological model of Cheng and co-workers [Cheng, E.; Cole, M. W.; Saam, W. F.; Treiner, J. Phys. ReV. Lett. 1991, 67, 1007. Cheng, E.; Cole, M. W.; Saam, W. F.; Treiner, J. Phys. ReV. B 1993, 48, 18214] to predict accurately the evolution of wetting temperature with substrate strength. In addition, we consider the extent to which select prewetting data from different wall potentials collapse when scaled by either wetting or prewetting critical point properties. I. Introduction Fluids routinely interact with solid substrates in industrial processes, natural phenomena, and throughout our daily experiences. Examples vary from the use of liquid detergents to clean stained fabric to the collection of rainwater onto tree leaves. The behavior of a fluid in the vicinity of a surface depends qualitatively upon the relative strengths and ranges of the fluid-fluid and fluid-substrate interactions and the structural characteristics of the substrate.1,2 It is the influence of the intermolecular interactions that we focus on here. Specifically, we use Monte Carlo simulation to investigate the evolution of prewetting phase behavior with substrate strength. Prewetting transitions are linked to first-order wetting points.3,4 The wetting point is typically associated with the location along the bulk liquid-vapor saturation line at which fluid adsorption switches from partial to complete wetting. At temperatures below the wetting point (partial wetting), fluids form droplets with nonzero contact angle, whereas above the wetting temperature (complete wetting), fluids spread evenly over a surface. The wetting point also represents a location on the phase diagram where one observes coexistence between a thin and infinitely thick fluid layer adsorbed to the substrate. Prewetting transitions represent an extension of this saturation point to state conditions off of the bulk vaporization line. It follows that prewetting saturation lines terminate at a wetting point. Phase equilibrium is characterized by coexistence between thin and thick (but finite) substrate films at a chemical potential (pressure) below the bulk saturation value at a given temperature. Akin to bulk vapor-liquid saturation lines, prewetting lines terminate at critical end points. These phase transitions are often observed with substrates of moderate strength. Experimentally, prewetting * To whom correspondence should be addressed. E-mail: jerring@ buffalo.edu.
transitions have been identified with systems such as 4He on cesium5,6 and mercury on sapphire.7 The goal of the current study is to develop a quantitative understanding of how prewetting saturation properties evolve with the relative strength of the substrate-fluid interaction. To probe this issue, we investigate the properties of a model in which a fluid with fixed characteristics contacts a substrate of variable strength. We examine the evolution of the two terminal points (wetting and prewetting critical temperatures) that define prewetting lines with variation of substrate strength and thereby identify the region within the temperature-wall strength space where prewetting coexistence is observed. We also investigate the ability of corresponding states approaches to describe prewetting behavior over a broad range of conditions. The work of Pandit et al. established general qualitative trends for the adsorption of fluids onto substrates of variable strength.8 They suggested that systems with substrate-fluid interactions of intermediate strength will exhibit wetting points, with the wetting temperature increasing with decreasing interaction strength. Their analysis also revealed that for strong enough substrates, the wetting temperature drops below a roughening temperature, signaling the onset of layering behavior. Conversely, as the substrate-fluid interaction is weakened, the wetting point approaches a critical wetting point, above which prewetting transitions cease to exist. Nakanishi and Fisher9 and later Ebner and Saam10 further detailed the surface phase behavior of fluids in the vicinity of a bulk critical point. These studies illustrated important differences in the expected behavior of systems with short- and long-ranged substrate-fluid interactions. A phenomenological model introduced by Cheng et al.11,12 led to a significant advance in the ability of scientists and engineers to quantitatively predict the location of wetting transitions. The model stems from a balance between the energy
10.1021/jp803458x CCC: $40.75 2008 American Chemical Society Published on Web 07/31/2008
12906 J. Phys. Chem. C, Vol. 112, No. 33, 2008 gain from an infinitely thick liquid layer in contact with a surface and the energy cost of constructing interfaces between the liquid slab and both the vapor and solid substrate. Implementation of the approach requires information regarding the substrate-fluid interaction potential and bulk saturation properties only. This simple model has proven effective in guiding experimentalists to state conditions where prewetting coexistence is likely to be found. For example, the location of the prewetting saturation line for 4He on cesium by Rutledge and Taborek5 was assisted by the Cheng et al. heuristic model. Gatica et al. recently reported the use of this approach to describe the wetting behavior of water on graphite.13 The validity of the model has been investigated via molecular simulation. Studies by Cole and co-workers14,15 and Ancilotto and Toigo16 have found the model to have semiquantitative accuracy, with the largest deviations from simulation observed for weak substrates. We further analyze the reliability of this model here. In this work, we examine the prewetting behavior of a Lennard-Jones system in which monatomic particles interact with a planar structureless substrate. We take the parameter set introduced by Ebner and Saam4 for the adsorption of argon on solid carbon dioxide as a reference system and investigate variation of prewetting properties upon a change in the depth of the attractive well within the substrate-fluid interaction potential. The Ar-CO2 reference system has been studied extensively via molecular simulation.17–23 The important work of Monson and co-workers established the existence of prewetting coexistence within this system.17–19While subsequent work helped refine the quantitative description of this system, accurate depiction of the prewetting saturation line near the wetting point proved difficult. We recently showed23 how one can overcome this difficulty through implementation of transition matrix Monte Carlo techniques.24–28 Specifically, we demonstrated how one can use a grand canonical version of this approach29,30 to locate prewetting coexistence at nearwetting temperatures and thereby determine the wetting temperature accurately. Determination of prewetting critical temperatures has also proven to be a challenging task. We address this issue here through examination of the boundary tension31 in the vicinity of the critical point. This quantity is analogous to the surface tension for bulk liquid-vapor transitions and, as such, vanishes at the prewetting critical point. As demonstrated earlier,32 grand canonical transition matrix Monte Carlo simulation can be combined with finite-size scaling techniques33 to compute the boundary tension. Adoption of this approach enables us to obtain system-size-independent values of the prewetting critical temperature. Collectively, the suite of Monte Carlo simulation techniques employed here provides us with a means to determine prewetting saturation properties in a quantitatively rigorous manner. The paper is organized as follows. In the following section, we describe the model examined in this work. We then introduce the simulation tools employed here and provide details related to our calculations. Next, we present our simulation results and discuss the ability of corresponding states approaches to characterize prewetting phase behavior. Finally, the salient conclusions are reviewed. II. Model System We examine the behavior of a collection of monatomic particles in the presence of a planar structureless substrate. The energy of interaction uff between any two fluid particles separated by a distance r is given by the truncated LennardJones 12-6 potential34
uff(r) )
{
Sellers and Errington
[( σr ) - ( σr ) ]
4ε
12
6
0
for r < rc for r > rc
(1)
where ε and σ are energy and size parameters, respectively, and rc ) 2.5σ is the cutoff distance. The energy of interaction uwf between an adsorbing wall and a fluid particle separated by a distance z is given by the Lennard-Jones 9-3 potential
uwf(z) )
[ ( ) ( )]
σw 2 σw 9 2π Fwσw3εw 3 15 z z
3
(2)
with Fwσw3 ) 0.988 and σw/σ ) 1.0962. We study the behavior of the system at multiple values of the substrate strength εw/ε. When εw/ε ) 1.2771, the parameter set matches that introduced by Ebner and Saam4 to describe the adsorption behavior of argon on solid carbon dioxide. This model system has been studied by a number of groups,17–23 and we use it here as a point of reference. More specifically, we describe the wall strength in terms of R ) εw/εw,Ar-CO2 ) (εw/ε)/1.2771. From this point forward, all quantities are made dimensionless using ε and σ as characteristic energy and length scales, respectively. For example, temperature and volume are reduced by ε/k and σ3, respectively. III. Molecular Simulation Approach We determine interfacial properties via analysis of the surface density dependence of the surface excess free energy at a given temperature, which is obtained using grand canonical transition matrix Monte Carlo simulation. Computation of a number of interfacial properties of interest requires additional information regarding the corresponding bulk system at the same temperature. In what follows, we briefly describe the general grand canonical transition matrix Monte Carlo simulation approach employed to obtain bulk and interfacial thermodynamic properties. A more detailed description of these methods can be found within the reports from our earlier studies.23,29,32,35 Grand Canonical Transition Matrix Monte Carlo. Grand canonical transition matrix Monte Carlo (GC-TMMC)29 is a general method for determining the density dependence of the free energy of a fluid at a given temperature. Simulations are conducted in a standard grand canonical ensemble36 where the volume V, chemical potential µo, and temperature T are held constant and the particle number N (density) and energy E fluctuate. During a Monte Carlo simulation, transitions between states of different density are attempted, and the probability of completing these transitions is tracked using a simple bookkeeping scheme.27,28 At regular intervals during a simulation, this information is used to obtain an estimate of the density probability distribution, which is subsequently used to bias the sampling37 to low probability densities. Over time, all densities of interest are sampled adequately. The end result is an efficient self-adaptive method for determining the density probability distribution over a density range of interest (e.g., a range that includes the densities of two potentially coexisting phases). The key quantity extracted from a GC-TMMC simulation is the particle number (density) probability distribution Π(N). From a thermodynamics viewpoint, this function provides the density dependence of the Helmholtz free energy at a given temperature.38,39 With such information, one can use standard thermodynamic relationships to determine thermophysical properties of interest. The density dependence of a system is often explored with the aid of histogram reweighting.40 This technique provides the particle number probability distribution at a chemical potential µ other than that originally simulated
Influence of Substrate Strength on Wetting Behavior
ln Π(N;µ) ) ln Π(N;µo) + Nβ(µ - µo)
J. Phys. Chem. C, Vol. 112, No. 33, 2008 12907
(3)
where β ) 1/kT (k is Boltzmann’s constant). We now describe how to analyze particle number probability distributions obtained from bulk and interfacial simulations to obtain bulk saturation and surface properties. Bulk Properties. To obtain bulk properties at a temperature of interest, we first employ GC-TMMC simulation to determine Π(N) over a range of particle numbers. Since we wish to extract both single-phase and liquid-vapor saturation properties from Π(N), we sample a range of particle numbers that includes those states characteristic of both the vapor and liquid phases. The bulk density Fb and pressure pb of a single-phase system at chemical potential µ and temperature T are given by29
Fb(µ, T) ) V-1
∑ NΠ(N;µ) N
[∑
pb(µ, T) ) TV-1 ln
N
Π(N;µ) Π(0;µ)
(4)
]
(5)
Note that when Π(N) is normalized, pb(µ,T) ) -TV-1 ln[Π(0;µ)]. The statistical mechanical formalism that provides the connection between the thermodynamic pressure and the density probability distribution can be found elsewhere.29,41 When the chemical potential is sufficiently close to a saturation value, Π(N) is bimodal, and the summations above are restricted to those regions of N that correspond to the phase of interest. The coexistence point is identified by using histogram reweighting to determine the chemical potential µbsat that results in equal pressures (psat b ) within the liquid and vapor phases. Part of the analysis below requires the bulk liquid-vapor surface tension. We used an expanded ensemble area sampling scheme42 to determine this quantity. Details regarding the simulations performed here are identical to those presented in ref 42. Several bulk saturation properties are collected in Table 1. Interfacial Properties. To calculate interfacial properties, we perform a GC-TMMC simulation using a rectangular parallelepiped box with periodic boundary conditions applied in the x and y directions. The system is closed at the two ends of the nonperiodic z direction by square adsorbing and hard walls of cross-sectional area A. Simulations are typically completed at the bulk fluid liquid-vapor saturation chemical potential µsat b for the temperature T of interest. Again, the key quantity extracted from such a simulation is the particle number probability distribution Π(N). The negative of the logarithm of the probability distribution provides the so-called surface excess free energy as a function of the surface density.43,44 The qualitative features that this curve adopts under different wetting scenarios have been described in a number of theoretical studies.43–46 The excess surface density Γ and interfacial tension γ at chemical potential µ and temperature T of a single-surface-phase fluid in contact with the adsorbing wall are given by35
Γ(µ, T) ) A-1
∑ NΠ(N;µ) - HFb(µ, T) - Γhw(µ, T) N
[∑
γ(µ, T) ) -TA-1 ln
N
(6)
]
Π(N;µ) + Hpb(µ, T) - γhw(µ, T) Π(0;µ) (7)
where Γhw and γhw represent the excess surface density and interfacial tension for the fluid in contact with a hard wall. These quantities can be extracted from simulations with hard walls placed at both ends of the simulation cell. In this case, Γ ) Γhw
TABLE 1: Bulk Saturation Properties T
µsat b
0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15
-4.1967(10) -4.0633(5) -3.9444(7) -3.8377(8) -3.7421(6) -3.6570(4) -3.5818(4) -3.5154(3) -3.4572(4) -3.4063(3) -3.3619(2) -3.3239(4)
3 F 3 psat b × 10 b,v × 10
0.557(1) 1.279(1) 2.588(4) 4.737(4) 8.020(8) 12.74(1) 19.21(1) 27.76(1) 38.75(4) 52.53(4) 69.59(1) 90.54(7)
0.938(2) 2.011(1) 3.831(5) 6.678(6) 10.88(1) 16.80(1) 24.93(1) 35.92(1) 50.79(7) 71.26(7) 101.3(1) 155.7(4)
Fb,l
γlv
∆hvap
0.8581(12) 0.8377(7) 0.8163(8) 0.7934(6) 0.7683(6) 0.7427(5) 0.7152(3) 0.6858(9) 0.6523(6) 0.6144(6) 0.5673(4) 0.4972(10)
0.988(16) 0.888(12) 0.781(19) 0.696(11) 0.595(9) 0.503(13) 0.408(2) 0.322(3) 0.243(7) 0.164(5) 0.098(2) 0.032(1)
6.30 6.13 5.94 5.72 5.48 5.20 4.86 4.48 3.99 3.37
and γ ) γhw, and the hard-wall quantities are obtained from eqs 6 and 7. Although we explicitly computed these quantities here, they did not contribute significantly to the values for Γ and γ reported below. Note that within the expressions above, we have taken the dividing surface to be the z ) 0 plane defined by eq 2. When the system is in the vicinity of a prewetting saturation point, Π(N) is bimodal, and the summations above are restricted to those regions of N that correspond to the phase of interest (thin or thick).23 Akin to the bulk case, the coexistence point is identified by using histogram reweighting to determine the sat chemical potential µpw that results in equal interfacial tensions sat (γ ) within the thin and thick surface phases. Boundary Tension. Prewetting transitions are inherently twodimensional in nature. As a result, the boundary between two coexisting phases is a one-dimensional structure (a line), defined by the edge-on meeting of a thin and a thick film. Such a boundary is associated with an inhomogeneity in density, which gives rise to an excess free energy. The magnitude of this quantity, per unit length of the interfacial line, defines the boundary tension τ. Conceptually, this property is analogous to the surface tension between two bulk phases.31 The boundary tension was used to locate prewetting critical points. More specifically, we identified the prewetting critical temperature Tpwc as the temperature for which τ(Tpwc) ) 0. The finite-size scaling formalism of Binder33 was used to determine this quantity. In this approach, a true boundary tension value τ at a given temperature is related to a set of apparent systemsize-dependent boundary tension values τL, obtained through a series of molecular simulations with varying substrate area, A ) L2, through the scaling relationship
1 ln L τL ) c1 + c2 +τ L L
(8)
where c1 and c2 are constants. The expression suggests that the term τL becomes linear in the scaling variable ln(L)/L as the system size approaches infinity. The formalism enables one to extrapolate the true infinite-system-size interfacial tension from a series of finite-system-size calculations.32 Apparent system-size-dependent boundary tensions are directly related to the free-energy barrier a system must traverse when transforming from one surface phase to another. For a prewetting transition, the relevant order parameter that connects two coexisting phases is the surface density. It follows that the magnitude of the free-energy barrier FL can be written in terms of the probability Π(N) of finding the system with a surface density corresponding to one of the two saturated phases (Πthin and Πthick), relative to the probability of the least likely macrostate along the path connecting the saturated states (Πmin)29,30,32
12908 J. Phys. Chem. C, Vol. 112, No. 33, 2008
1 βFL ) (ln Πthin + ln Πthick) - ln Πmin 2
Sellers and Errington
(9)
The system-size-dependent boundary tension is given by the magnitude of the free-energy barrier per unit length of the interfacial line, τL ) FL/2L, where the 2 appears in the denominator due to the presence of two interfaces within a given simulation cell. Quantification of the prewetting critical temperature requires one to enumerate τ over a range of temperatures near Tpwc. This information is obtained by first collecting a joint particle number and energy probability distribution Π(N,E) during a GC-TMMC sat simulation conducted at estimates for Tpwc and µpw (Tpwc). We use a transition matrix approach to capture the density dependence and the visited states technique to obtain the energy dependence of this distribution. A detailed description of this approach is provided in ref 29. Once Π(N,E) has been constructed at the aforementioned temperature and chemical potential values, histogram reweighting40 is used to obtain the corresponding distribution at some nearby temperature T (and alternate chemical potential if desired). This joint distribution is then collapsed to a particle number probability distribution Π(N) specific to the nearby temperature T. Evaluation of the free-energy barriers FL and apparent boundary tensions τL at this temperature then follow from Π(N) as described above.
Figure 1. Representative particle number probability distributions for the model system studied in this work. (a) Distributions for various substrate strengths at T ) 0.75. Curves from bottom to top span from R ) 0.7 to 1.3 in increments of 0.1. (b) Distributions for various temperatures with R ) 0.90. Curves from bottom to top span from T ) 0.65 to 1.00 in increments of 0.05.
IV. Simulation Details Bulk simulations of the truncated Lennard-Jones fluid were performed over a temperature range spanning from T ) 0.600 to 1.150 in increments of 0.025. GC-TMMC simulations were conducted in a cubic cell with a volume of V ) 512. Other aspects of the simulations are similar to those described in ref 29. We examined substrate strengths spanning from R ) 0.5 to 1.5. The majority of interfacial simulations were performed with a box of area A ) 81 and height that spanned between H ) 20 and 80. The maximum particle number used in the GC-TMMC simulations varied between 800 and 3000. Simulation boxes with areas of A ) 144, 256, and 400 were used to compute boundary tension values. The Monte Carlo move mix included standard particle displacements and additions/deletions as well as aggregation-volume-bias47–49 displacements (AVBMC2) and additions/deletions. Statistical uncertainties were determined by performing four independent sets of simulations. The standard deviation of the results from the four simulation sets is taken as an estimate of the statistical uncertainty. V. Results and Discussion Molecular Simulation Results. The primary objective of this study was to characterize the influence of temperature and substrate strength on the prewetting phase behavior exhibited by the model system outlined above. To accomplish this task, we performed a series of molecular simulations that enabled us to trace the loci of wetting and prewetting critical temperatures (Tw and Tpwc) as a function of substrate strength. To describe the phase behavior at a given set of conditions, we calculated and analyzed the particle number probability distribution, which provides the surface excess free energy as a function of surface density.35 The relationship between this curve and several macroscopic wetting properties (e.g., contact angle, interfacial tension, prewetting saturation densities) has been described by a number of authors.43–46
Figure 1a and b provides particle number probability distributions for several wall strengths at a temperature of T ) 0.75 and for various temperatures at R ) 0.9, respectively. Each of the distributions in Figure 1 is evaluated at µ ) µsat b , and therefore, in the limit of large N, wherein the influence of the substrate on a particle diminishes, the curves tend toward a limiting “plateau” value. One can use these curves to map the system behavior to one of three general categories: partial wetting conditions (T < Tw), prewetting conditions (Tw < T < Tpwc), and complete wetting (and super prewetting critical) conditions (T > Tpwc). Partial wetting conditions correspond to a situation in which the plateau probability is smaller than that at low surface densities. In this case, Π(N) can be used to extract the spreading coefficient.35 Prewetting scenarios are described by a plateau probability that is larger than that at low surface densities along with a change in the curvature of Π(N) at intermediate N. Note that this condition prevails when the system is in the vicinity of a first-order wetting transition only, which is characterized by a free-energy barrier between low-density (thin) and high-density (thick) surface phases. Finally, when T > Tpwc, Π(N) monotonically increases to its limiting value with no change in curvature. The curves in Figure 1a indicate that the wetting behavior evolves from partial (T < Tw) to complete (T > Tpwc) wetting conditions as the substrate strength is increased at constant temperature. The data provided in Figure 1b suggest an analogous evolution of wetting behavior with increasing temperature at constant substrate strength. To quantitatively determine the wetting temperature for a given wall strength, we examined the relationship between the sat chemical potential difference ∆µ ) µsat pw - µb and temperature. These data are presented in Figure 2 for various wall strengths. Theoretical arguments50,51 suggest that ∆µ ∝ (T - Tw)3/2 along the prewetting saturation line. Wetting temperatures were deduced via a fit to this functional relationship using the four points closest to Tw and are collected in Table 2. As is expected, the wetting temperature increases with decreasing substrate strength. At wall strengths stronger than R ) 1.0, layering/
Influence of Substrate Strength on Wetting Behavior
Figure 2. Temperature dependence of the chemical potential difference sat ∆µ ) µpw - µsat b . Squares, down triangles, ×’s, circles, pluses, up triangles, and diamonds correspond to substrate strengths of R ) 0.80, 0.90, 0.95, 1.00, 1.05, 1.10, and 1.20, respectively. Filled symbols signify the prewetting critical point. Broken curves represent extrapolations to the wetting point (see text for details).
TABLE 2: Wetting and Prewetting Critical Temperatures R
Tw
Tpwc
0.7 0.8 0.9 1.0 1.1 1.2 1.3
0.958(6) 0.846(3) 0.727(6) 0.582(3)
1.070(3) 1.006(2) 0.944(3) 0.875(2) 0.805(2) 0.731(3) 0.660(4)
crystallization behavior was observed at temperatures above the extrapolated wetting point. Under these conditions, the probability distribution showed patterns different from those described above (e.g., oscillatory patterns) and proved more difficult to converge than that with an amorphous fluid in contact with the substrate. Although these findings pointed to the emergence of interesting behavior at these conditions, we did not study this region of the phase diagram in detail. The description provided here regarding layering/crystallization is purely qualitative in nature. Quantitative prediction of properties within this region would require attention to additional simulation details (e.g., allowing the wall lengths in the x and y directions to fluctuate so that crystalline layers can fully relax). We also identified the high-temperature terminal point along various prewetting saturation lines. This task was accomplished by finding the first temperature Tpwc at which the boundary tension τ evaluated to zero (from below). An illustrative example of the information examined to determine the prewetting critical temperature for the R ) 0.9 case is presented in Figure 3. For each substrate strength of interest, GC-TMMC simulations were performed with three system sizes (A ) 144, 256, 400) at an estimate for Tpwc (T ) 0.940 for the R ) 0.9 case). The particle number probability distributions from these simulations (see Figure 3a) were used to obtain apparent boundary tensions τL for each of the system sizes. Finite-size scaling (see Figure 3b) was subsequently employed to extrapolate to the true (infinitesystem) boundary tension τ at the estimated prewetting critical temperature. As noted within the Molecular Simulation Ap-
J. Phys. Chem. C, Vol. 112, No. 33, 2008 12909
Figure 3. (a) Particle number probability distributions for the R ) 0.9 system at the near-critical temperature of T ) 0.94. Curves from bottom to top represent system sizes of A ) 144, 256, and 400. (b) Finite-size scaling of the apparent boundary tension. Squares, circles, diamonds, and left triangles correspond to temperatures of T ) 0.935, 0.940, 0.945, and 0.950, respectively. Open symbols represent finitesize data, and closed symbols represent extrapolations to the infinitesystem limit. (c) Infinite-system boundary tension as a function of temperature.
proach section above, density-specific energy distributions were collected during each GC-TMMC run. This information enabled us to construct particle number probability distributions over a range of temperatures in the vicinity of that originally simulated and, therefore, predict the temperature dependence of the boundary tension at conditions close to Tpwc. Figure 3c provides the τ(T) data obtained for the R ) 0.9 case. Adopting this approach enabled us to determine infinite-system-size prewetting critical temperatures. Numerical values for Tpwc at several R’s are collected in Table 2. Once a value for Tpwc was obtained, we extrapolated subcritical data for ∆µ(T) and γsat(T) to the prewetting critical point to estimate ∆µpwc and γpwc, respectively. The prewetting critical excess surface density Γpwc was approximated from the point at which Π(N) for the A ) 400 system size exhibited a minimum within the interfacial region. Figure 4 provides a summary of prewetting phase behavior for the model Lennard-Jones system studied in this work. Our results indicate that both the wetting and prewetting critical temperatures increase with decreasing substrate strength, with a sharper increase in the wetting temperature. It follows that the length of the prewetting saturation line (defined by the temperature difference ∆T ) Tpwc - Tw) decreases with decreasing wall strength. Mean-field theoretical work by Ebner and Saam suggests that the wetting transition should remain first-order upon approach to the bulk critical temperature Tc when the substrate-fluid interaction potential has a sufficiently long range.10 Given that the surface potential considered here includes a long-range van der Waals attraction, we should expect the wetting transition to remain first-order to Tc. We found that both the wetting and prewetting critical lines were wellrepresented by the scaling form (T - Tc) ∝ R3/2. The continuous lines within Figure 4 result from curve fits with this functional
12910 J. Phys. Chem. C, Vol. 112, No. 33, 2008
Figure 4. Surface phase diagram for the model Lennard-Jones system studied here. Diamonds and circles represent wetting and prewetting critical temperatures, respectively. The broken curves stem from a fit to a scaling form detailed within the text. Statistical uncertainties are smaller than the size of the symbol.
Figure 5. Temperature dependence of saturated excess surface densities. Symbols are equivalent to those used in Figure 2.
relationship. The two curves meet at the bulk critical temperature52 of Tc ) 1.1876 and wall strength of R ≈ 0.47, where we expect the prewetting region to vanish. Particle number probability distributions (free-energy curves) provide additional support for this conclusion. The Π(N) curve for T ) 1.10 and R ) 0.5 clearly indicates a partial wetting scenario with firstorder wetting behavior. Analogous data collected at T ) 1.15 do not provide the same level of clarity. Although the behavior at this temperature remains first-order, we could not discern whether the system resided within the partial wetting regime. The convergence of Π(N) to its limiting “plateau” value was slow and proved difficult to extrapolate with confidence. Overall, the simulation data collected here suggest that prewetting phase transitions are observed within a triangular-shaped region within the temperature-wall strength plane, enclosed by wetting, critical prewetting, and layering/crystallization boundaries. The relatively small size of this region, coupled with the small values documented for the chemical potential differences ∆µ, contributes to the difficulties often encountered in locating prewetting transitions (via both experiment and simulation).2 Figure 5 shows saturated excess surface densities for the thin and thick phases. For temperatures close to the wetting point, the density of the thick phase grows rapidly with decreasing temperature. In the absence of a gravitational field, the high-
Sellers and Errington
Figure 6. Temperature dependence of the saturated interfacial tension (main panel) and enthalpy of thinning (inset). Symbols are equivalent to those used in Figure 2.
density layer becomes infinitely thick at the wetting temperature. The rapidly increasing density of the thick phase renders the location of prewetting saturation intractable via simulation at temperatures close to Tw. This complication becomes progressively more severe as one moves to lower values of R, wherein the temperature gap between the wetting and prewetting critical point decreases and the wetting temperature increases to higher values, leading to enhanced density fluctuations (broader peaks) near Tw. Figure 6 provides prewetting saturation values for the interfacial tension. Our results for the R ) 1.0 case are in reasonably good agreement with those from previous studies. Fan and Monson19 and Sweatman22 provided T ) 0.88 estimates of γsat ) -0.180 and -0.165, respectively. Although our results indicate that the T ) 0.88 isotherm lies slightly above the prewetting critical point (Tpwc ) 0.875), one can modestly extrapolate the saturation line to obtain an estimate for comparison of γsat ) -0.170 at this temperature. Analysis of the interfacial tensions indicates that ln(-γsat) is nearly linear in T-1. The temperature dependence of the saturated interfacial tension is related to the molar “heat of thinning” ∆hpw ) hthin - hthick, which we estimated as follows
( )(
∆hpw ) T∆spw ) -T
dγsat -1 -1 Γthin - Γthick ) dT
(10)
where ∆spw is the molar entropy of thinning. The derivative in the expression above was approximated using a center-point difference with (dγsat/dT) ≈ -γsatT-2[∆ln(-γsat)/∆(1/T)]. Values for ∆hpw are displayed within the inset to Figure 6. As one might intuitively expect, the magnitude of this quantity is smaller than the heat of vaporization ∆hvap for the bulk fluid at the same temperature (see Table 1). Before ending the discussion regarding our molecular simulation results, we comment briefly on the influence of system size on prewetting saturation properties. As is apparent from the discussion above, the only property for which we rigorously accounted for finite-size effects was the critical prewetting temperature. Preliminary calculations reveal statistically significant variation in the saturation properties with system size. For example, an extrapolation of µsat pw to the infinite system limit shows an increase in this quantity of 0.002 relative to the A ) 81 system size for conditions of R ) 1.0 and T ) 0.70. The magnitude of this change is approximately a factor of 10 larger
Influence of Substrate Strength on Wetting Behavior
J. Phys. Chem. C, Vol. 112, No. 33, 2008 12911
Figure 7. Comparison of simulation (diamonds) and CCST (squares) estimates for the dependence of substrate strength on the wetting temperature. Simulation data are equivalent to those displayed in Figure 4. sat than the statistical uncertainty in µsat pw. This increase in µpw leads to a 5% increase in the value of ∆µ, which, provided the shift in ∆µ remains positive at all T, suggests that the infinite system wetting temperature should be slightly higher than the value reported above. Bulk properties away from Tc exhibit much less dependence on system size. This difference can be understood by recognizing that saturation properties generally scale with system size as L-d, where d is the dimensionality of the phase transition. We plan to address the role of finite-size effects in greater detail in future work. Corresponding States Analysis. A number of ideas have been put forth to describe wetting phase behavior in a universal manner. Cheng et al. (CCST) have introduced a simple phenomenological model for predicting the wetting temperature.11,12 Recall that the wetting transition is characterized by coexistence between a thin surface film and an infinitely thick film at a point along the bulk liquid-vapor saturation line (Tw, sat µsat b ) µpw). Within the CCST model, the wetting temperature stems from a balance between the potential energy gain from the interaction of an infinitely thick fluid slab with a solid substrate and the energy cost of forming solid-liquid and liquid-vapor interfaces. The wetting temperature is given by the implicit expression
[Fb,l(Tw) - Fb,v(Tw)] ∫zmin uwf(z)dz ) -2γlv(Tw) (11)
[
∞
]
where zmin represents the position at which uwf reaches a minimum value, γlv is the liquid-vapor surface tension, and Fb,l and Fb,v are the saturated liquid and vapor densities of the bulk fluid, respectively. Provided one has the ability to deduce the substrate-fluid interaction potential uwf and bulk saturation properties of the fluid (via experiment, simulation, or theory), eq 11 provides a general means for estimating Tw. Gatica et al. recently detailed the approximations that accompany this approach.13 Here, we examine the extent to which our simulation-based values for Tw agree with those from the CCST model. Inclusion of the substrate-fluid potential used within this study (eq 2) leads to the following simplified version of eq 11
R(Tw) ) 1.110γlv(Tw)/[Fb,l(Tw) - Fb,v(Tw)]
(12)
Data for the bulk properties are included in Table 1. Figure 7 contains CCST and simulation estimates for the evolution of the wetting temperature with substrate strength. The CCST model produces a much weaker rate of change in wetting
Figure 8. Scaled prewetting phase diagram. Symbols are equivalent to those used in Figure 2. Statistical uncertainties are smaller than the size of the symbol.
temperature with substrate strength relative to simulation. The CCST model indicates that the first-order wetting line extends to low substrate strength with Tw f Tc as R f 0. This picture is in contradiction to that depicted by Ebner and Saam10 in which the wetting line terminates at nonzero R, with a locus of firstorder drying transitions emerging at low substrate strength. Nonetheless, the CCST model provides a reasonable prediction for Tw at intermediate values of R. As shown earlier by Curtarolo et al.14 and Ancilotto and Toigo,16 CCST predictions underestimate Tw for weak substrates and, conversely, overestimate Tw for strong substrates. Bonn and Ross have demonstrated that experimental prewetting data for ∆µ(T) from a seemingly diverse set of systems roughly collapse onto a single curve when ∆µ and T are appropriately scaled by Tw.2 More specifically, they found that
-
( )
∆µ(T) 1 T - Tw = Tw 2 Tw
3/2
(13)
which stems from modification of the expression used to deduce Tw above (see Figure 2). Our data for several R values are collected in Figure 8. Remarkably, all of the data collapse onto the corresponding states curve presented by Bonn and Ross. The fact that each of the data sets shares the same power law scaling exponent of 3/2 is not surprising as this value is expected due to the z-3 van der Waals term present in the substrate-fluid interaction potential.50,51 However, the observation that all substrate strengths share the same proportionality constant (1/2 in eq 13) is perhaps not generally expected. The experimental data presented by Bonn and Ross displayed greater variability in this value. Consideration of a broader range of fluid systems may help in understanding the extent to which this constant is universal. Bonn and Ross also noted that the prewetting critical points for several systems were found at -∆µ/Tw ≈ 0.03. We do not arrive at the same conclusion. Our values for -∆µ/Tw at Tpwc span from approximately 0.03 to 0.2. In fact, taking -∆µ/Tw to be a constant value at Tpwc implies that the quantity (Tpwc - Tw)/Tw would remain invariant with substrate strength (recall eq 13). This result suggests that ∆T ) Tpwc - Tw would increase with increasing Tw (decreasing R), which is in contradiction to the picture presented in Figure 4. We also examined the possibility of scaling prewetting saturation data by the upper prewetting critical end point, akin to the scaling performed routinely for bulk fluids. Figure 9
12912 J. Phys. Chem. C, Vol. 112, No. 33, 2008
Sellers and Errington these properties. Collectively, these corresponding states results provide a basis to predict prewetting coexistence properties with either wetting or prewetting critical point information alone. Finally, our results for the latent heat of the prewetting transition suggest that this quantity evaluates to approximately 60% of the bulk enthalpy of vaporization at the same reduced temperature. The results presented here point to several directions for future research. Much of the corresponding states analysis presented above would benefit from a broader range of well-characterized systems to consider. Inclusion of systems with greater diversity in the form of both the fluid and substrate interaction potential would be valuable. The methods employed here proved robust and hold promise for constructing phase diagrams of interest in a rapid and efficient manner. Finally, a more detailed analysis of the influence of finite-size effects on prewetting saturation properties would be helpful.
Figure 9. Corresponding states prewetting phase diagram. Symbols are equivalent to those used in Figure 2. Broken curves provide equivalent bulk properties along the vaporization line. Statistical uncertainties are smaller than the size of the symbol.
contains reduced “thinning” curves for several values of R. We find that the data approximately collapse onto a single curve, with the magnitude of the slope increasing slightly with increasing R. Interestingly, we also find that the prewetting saturation data track closely the reduced bulk vaporization line sat characterized by the relationship between psat b /pb,c and T/Tc. The inset to Figure 9 provides the heat of thinning as a function of reduced temperature. Again, we find that all of the prewetting data collapse onto a single curve. Comparison to the bulk equivalent reveals that ∆hpw ≈ 0.6∆hvap at the same reduced temperature. VI. Conclusions We have examined the influence of substrate strength on prewetting phase behavior. Grand canonical transition matrix Monte Carlo simulation was used to determine saturation properties of interest, including excess surface densities, interfacial tensions, chemical potentials, and latent enthalpies of thinning. Wetting temperatures were obtained from consideration of the difference in bulk and prewetting saturation chemical potentials along the prewetting line. The temperature dependence of the boundary tension was examined to deduce prewetting critical temperatures. Overall, our results indicate that prewetting phase behavior resides within a triangular-shaped region within the temperature-wall strength plane. One vertex of the triangle is defined by intersection of the wetting and prewetting critical lines, which occurs at the bulk critical temperature and nonzero wall strength. The other two vertices are associated with the first temperature at which layering/crystallization occurs when moving along the wetting and prewetting critical lines. Analysis of the heuristic model of Cheng et al.11,12 indicated that this approach leads to accurate values for the wetting temperature at intermediate substrate strength only. For the Lennard-Jones system studied here, the model underestimated (overestimated) the wetting temperature for weak (strong) substrates. Following the work of Bonn and Ross,2 we examined the temperature dependence of the difference between bulk and prewetting saturation chemical potentials along several saturation lines. We found that these data collapse onto a single curve when appropriately reduced by the wetting temperature. We also found that interfacial tension-temperature data collapse onto a common curve when reduced by prewetting critical values for
Acknowledgment. We gratefully acknowledge the financial support of the Donors of the American Chemical Society Petroleum Research Fund (Grant No. 43452-AC5) and the National Science Foundation (Grant No. CTS-028772). Computational resources were provided in part by the University at Buffalo Center for Computational Research. References and Notes (1) Gelb, L. D.; Gubbins, K. E.; Radhakrishnan, R.; SliwainskaBartkowiak, M. Rep. Prog. Phys. 1999, 62, 1573. (2) Bonn, D.; Ross, D. Rep. Prog. Phys. 2001, 64, 1085. (3) Cahn, J. W. J. Chem. Phys. 1977, 66, 3667. (4) Ebner, C.; Saam, W. F. Phys. ReV. Lett. 1977, 38, 1486. (5) Rutledge, J. E.; Taborek, P. Phys. ReV. Lett. 1992, 69, 937. (6) Ketola, K. S.; Wang, S.; Hallock, R. B. Phys. ReV. Lett. 1992, 68, 201. (7) Yao, M.; Hensel, F. J. Phys.: Condens. Matter 1996, 8, 9547. (8) Pandit, R.; Schick, M.; Wortis, M. Phys. ReV. B 1982, 26, 5112. (9) Nakanishi, H.; Fisher, M. E. Phys. ReV. Lett. 1982, 49, 1565. (10) Ebner, C.; Saam, W. F. Phys. ReV. Lett. 1987, 58, 587. (11) Cheng, E.; Cole, M. W.; Saam, W. F.; Treiner, J. Phys. ReV. Lett. 1991, 67, 1007. (12) Cheng, E.; Cole, M. W.; Saam, W. F.; Treiner, J. Phys. ReV. B 1993, 48, 18214. (13) Gatica, S. M.; Johnson, J. K.; Zhao, X. C.; Cole, M. W. J. Phys. Chem. B 2004, 108, 11704. (14) Curtarolo, S.; Stan, G.; Bojan, M. J.; Cole, M. W.; Steele, W. A. Phys. ReV. E 2000, 61, 1670. (15) Bojan, M. J.; Stan, G.; Curtarolo, S.; Steele, W. A.; Cole, M. W. Phys. ReV. E 1999, 59, 864. (16) Ancilotto, F.; Toigo, F. Phys. ReV. B 1999, 60, 9019. (17) Finn, J. E.; Monson, P. A. Phys. ReV. A 1989, 39, 6402. (18) Finn, J. E.; Monson, P. A. Langmuir 1989, 5, 639. (19) Fan, Y.; Monson, P. A. J. Chem. Phys. 1993, 99, 6897. (20) Skowłowski, S.; Fischer, J. Phys. ReV. A 1990, 41, 6866. (21) Shi, W.; Zhao, X.; Johnson, J. K. Mol. Phys. 2002, 100, 2139. (22) Sweatman, M. B. Phys. ReV. E 2001, 65, 011102. (23) Errington, J. R. Langmuir 2004, 20, 3798. (24) Smith, G. R.; Bruce, A. D. J. Phys. A 1995, 28, 6623. (25) Wang, J. S.; Tay, T. K.; Swendsen, R. H. Phys. ReV. Lett. 1999, 82, 476. (26) Wang, J. S.; Swendsen, R. H. J. Stat. Phys. 2002, 106, 245. (27) Fitzgerald, M.; Picard, R. R.; Silver, R. N. Europhys. Lett. 1999, 46, 282. (28) Fitzgerald, M.; Picard, R. R.; Silver, R. N. J. Stat. Phys. 2000, 98, 321. (29) Errington, J. R. J. Chem. Phys. 2003, 118, 9915. (30) Errington, J. R. Phys. ReV. E 2003, 67, 012102. (31) Widom, B. Mol. Phys. 1999, 96, 1019. (32) Errington, J. R.; Wilbert, D. W. Phys. ReV. Lett. 2005, 95, 226107. (33) Binder, K. Phys. ReV. A 1982, 25, 1699. (34) (a) Lennard-Jones, J. E. Proc. R. Soc. London, Ser. A 1924, 106, 441. (b) Lennard-Jones, J. E. Proc. R. Soc. London, Ser. A 1924, 106, 463. (35) Grzelak, E. M.; Errington, J. R. J. Chem. Phys. 2008, 128, 014710. (36) Norman, G. E.; Filinov, V. S. High Temp. 1969, 7, 216. (37) Berg, B. A.; Neuhaus, T. Phys. ReV. Lett. 1992, 68, 9.
Influence of Substrate Strength on Wetting Behavior (38) Mittal, J.; Errington, J. R.; Truskett, T. M. Phys. ReV. Lett. 2006, 96, 177804. (39) Errington, J. R.; Truskett, T. M.; Mittal, J. J. Chem. Phys. 2006, 125, 244502. (40) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635. (41) Errington, J. R.; Shen, V. K. J. Chem. Phys. 2005, 123, 164103. (42) Errington, J. R.; Kofke, D. A. J. Chem. Phys. 2007, 127, 174709. (43) Indekeu, J. O. Int. J. Mod. Phys. B 1994, 8, 309. (44) Indekeu, J. O. Physica A 1992, 183, 439. (45) De Gennes, P. G. ReV. Mod. Phys. 1985, 57, 827.
J. Phys. Chem. C, Vol. 112, No. 33, 2008 12913 (46) Getta, T.; Dietrich, S. Phys. ReV. E 1998, 57, 655. (47) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8725. (48) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 11275. (49) Chen, B.; Siepmann, J. I.; Oh, K. J.; Klein, M. L. J. Chem. Phys. 2001, 115, 10903. (50) Hauge, E. H.; Schick, M. Phys. ReV. B 1983, 27, 4288. (51) Sullivan, D. E.; Telo da Gama, M. M. In Fluid Interfacial Phenomenam; Croxton, C. A., Ed.; Wiley: New York, 1986. (52) Wilding, N. B. Phys. ReV. E 1995, 52, 602.
JP803458X
