Molecular Dynamics Simulation of Heterogeneous Nucleation and

Pavel V. Ovcharov , Nikita P. Kryuchkov , Kirill I. Zaytsev , and Stanislav O. Yurchenko. The Journal of Physical Chemistry C 2017 121 (48), 26860-268...
0 downloads 0 Views 496KB Size
Download PDF
15784

J. Phys. Chem. C 2007, 111, 15784-15791

Molecular Dynamics Simulation of Heterogeneous Nucleation and Growth of Argon at Polyethylene Films† Roberto Rozas and Thomas Kraska* Physical Chemistry, UniVersity of Cologne, Luxemburger Strasse 116, D-50939 Ko¨ln, Germany ReceiVed: May 15, 2007; In Final Form: August 13, 2007

Heterogeneous nucleation and growth of supersaturated argon vapor at polyethylene surfaces is investigated using molecular dynamics simulation. The specific system is chosen as model for high wettable systems. The simulations are conducted in a nonequilibrium ensemble which includes heat transfer during the condensation process. Temporary density and temperature gradients are developed in the vapor phase. The gradual transition from nearly adsorption behavior close to the binodal to heterogeneous nucleation of a metastable vapor is investigated by stepwise varying the initial saturation of the vapor. By increasing the supersaturation along an isotherm a continuous transition from the layer-by-layer growth to the islands-on-layers one is observed. The results are compared to classical nucleation theory for the heterogeneous two- and three-dimensional growth models. We find for this system that a two-dimensional version of the classical heterogeneous nucleation theory is most suitable to describe the nucleation rate data vs saturation obtained from simulation over a wide range of supersaturation.

1. Introduction The understanding of heterogeneous nucleation is important for the development of materials1 as well as for the description of natural phenomena such as the formation of atmospheric aerosols.2 There are several theoretical attempts to model homogeneous and heterogeneous nucleation ranging from analytical approaches such as classical nucleation theory3-5 to thermodynamic density functional theory6,7 and molecular simulations including molecular dynamics8-12 and Monte Carlo methods.13-15 Classical nucleation theory provides an interpretation of experimental and simulation results by establishing a link between the relevant parameters of the phase transition such as the nucleation rate, formation energy, and size of the nucleus, and the macroscopic properties of the system such as density, surface energy, and contact angle. On the other hand, molecular simulations provide a direct method to access the properties of the studied system and further information about the growth and aggregation phenomena at the microscopic level. In the case of heterogeneous nucleation, there is an additional degree of freedom, which is the interaction strength between the substrate and the nucleating vapor. This interaction can vary over a wide range depending on the combination of the substrate material and the vapor, which leads to different nucleation and growth mechanisms. Zapadinsky et al.15 investigated the effect of a substrate on the energy and the cluster size distribution of the system for different cluster sizes and cluster substrate distances by canonical Monte Carlo simulation. Kholmurodov et al.12 studied the condensation of a vapor confined in a slit pore by means of molecular dynamics simulations of argon between atomic walls separated by a short distance. This is a special case of heterogeneous nucleation being different from the nucleation at a free surface.16 Kholmurodov et al. used an inert gas thermostat by applying Berendsen and Nose´-Hoover thermostats * Corresponding author. [email protected]. † Part of the “Keith E. Gubbins Festschrift”.

to half of the vapor atoms and maintaining the substrate temperature using a combination of a tethering force and a constraint mechanism. One interesting result of their work is that the heterogeneous nucleation rate at an atomic wall can be about 1 order of magnitude larger than that of an idealized continuous wall. Kimura and Maruyama17 investigated argon condensation on atomic surfaces. Their simulation box consists of an atomic surface on one side and a reflection wall on the other side in a distance of about 14 nm. This distance is actually rather close to the width of the slit pore of around 12 nm in the simulation of Kholmurodov et al.12 The nucleation rates obtained by Kimura and Maruyama agreed with the heterogeneous threedimensional growth model (eq 3) of the classical nucleation theory.17 Depending on the relative strength of the cohesive interaction between particles in the vapor phase and the adhesive interaction between particles and the substrate different growth mechanisms are possible. Clusters can be homogeneously formed in a vapor; they can growth on a surface as three-dimensional islands or as two-dimensional layers. In the classical heterogeneous nucleation theory the growth mechanism is directly related to the wettability of a surface, which can be measured in terms of the contact angle for cap-shaped clusters. High wettability conditions favor the formation of two-dimensional layers rather than threedimensional islands. The effect of the Lennard-Jones potential parameters for the attraction in the fluid and between the liquid and the solid surface on the wetting behavior has been studied by molecular simulations of a droplet equilibrated on an atomic crystalline surface.18 These simulation results show that the macroscopic Young’s equation is still valid for small droplets and that the cosine of the contact angle linearly depends on the ratio surf/, where surf denotes the depth of the integrated potential between the surface and adatom and  the depth of the potential between the ad-atoms. In a preceding work, we have investigated the deposition of platinum on polyethylene as model for metal-polymer sys-

10.1021/jp073713d CCC: $37.00 © 2007 American Chemical Society Published on Web 10/10/2007

Heterogenous Nucleation

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15785

tems.19 In that system, the interaction between the metal atoms is significantly larger than the interaction between the metal and the polymer groups. As a result, we observed the formation of three-dimensional metal clusters in the polymer surface. Here we investigate a system which, according to the study of Maruyama et al.,18 possesses a high wettability condition. Simulations at different vapor-phase densities and two different substrate temperatures are performed. The chosen state points cover a wide range of the metastable phase diagram of pure argon vapor. We use nonequilibrium molecular dynamics simulations in order to mimic the experiment rather than an idealized ensemble. Specially, the growth mechanism of clusters formed by condensation of supersaturated vapor on a structured substrate is investigated. 2. Method 2.1. Molecular Dynamics Simulation. The molecular dynamics (MD) computer programs are developed here for this specific investigation. All programs are verified by comparison with literature results for the used argon potential and the polyethylene model. For polyethylene we employ the model of Sumpter et al.20 This is a united-atom model where each site represents a CH2 group or a CH3 end group. The interaction between the sites in different chains and between sites in the same chain separated by at least three valence bonds is described by the Lennard-Jones potential. Three different internal degrees of freedom are included in the potential model: the covalent bond between adjacent sites is modeled by a Morse potential, bending motions are described by a three-body cosine-angle potential, and a four-body potential is used for internal torsional vibrations and conformational isomerism. Here, the Morse potential is approximated by a spring potential for more efficient calculations. All parameters for the polyethylene model are taken from the literature.20 As cutoff radius, the value 1.5σPE is used for the polyethylene sites, where σ is the Lennard-Jones diameter of a site. For argon and the cross interactions, a cutoff of 2.5σAr is used. The Lennard-Jones parameters for argon are σAr ) 0.3405 nm and Ar/kB ) 117.7 K, and for the polyethylene sites, they are σPE ) 0.39 nm and PE/kB ) 59.378 K. The cross interaction parameters are approximated by the LorentzBerthelot rules. The polymer film consists of 374 chains with 70 sites each, and the vapor phase contains 5000 argon atoms. To generate the initial configuration, the polyethylene film and the argon vapor are equilibrated separately. The polyethylene substrate is prepared at TPE ) 60 K and TPE ) 80 K in two steps: first, starting from a crystalline orthorhombic structure, the polymer is simulated 1 ns in an anisotropic NpT ensemble maintaining the diagonal components of the stress tensor at a value of 0. Then the simulation box is extended in the direction parallel to the orientation of the polymer chains to produce a film followed by another equilibration run of the polyethylene film in vacuum. Argon vapor is equilibrated at 200 K, which is above its critical temperature. Finally, the equilibrated polymer film at either 60 or 80 K is merged with the argon vapor phase at 200 K. In the proceeding simulation periodic boundary conditions are applied in all three dimensions but no barostat is applied, i.e., the box size is maintained. A Nose´-Hoover thermostat is applied to the polyethylene film only. In contrast to the work of Kholmurodov et al.12 here the temperature of the vapor phase is only influenced by interaction with the polymer film. Heat transfer takes place via the argon atoms, which cool down by collisions with the polyethylene film. This assembly aims to represent a corresponding experiment. The polymer film is rather thin, however, the application

Figure 1. Low-density part of the phase diagram of argon. Ar exp: experimental data of saturated vapor densities;24 LJ-EOS: saturated vapor densities calculated from a Lennard-Jones equation of state;25 LJ-2.5σ: saturated vapor densities calculated from the Lennard-Jones potential with a cutoff of 2.5σ; spinodal LJ-EOS: stability limit estimated from the Lennard-Jones equation of state. The crosses mark the initial state points of the vapor in the nucleation simulations.

of a thermostat to the film mimics a large heat reservoir similar to a surface of a macroscopic polymer block. The low surface temperature induces a supersaturation on the vapor near the surface and adsorbs the condensation heat. Simulations at six different argon vapor densities are performed (10, 20, 30, 40, 50, and 60 g/dm3 for TPE ) 80 K and 5, 10, 15, 20, 25, and 30 g/dm3 for TPE ) 60 K) where TPE denotes the temperature of the polymer substrate. These temperatures are lower than the glass transition temperature of the polyethylene film and also lower than the critical temperature of argon. The obtained lengths of the simulation box in the transversal directions are Lx ) 7.96 nm and Ly ) 7.95 nm for 80 K and Lx ) 7.93 nm and Ly ) 7.93 nm for 60 K. These values are a result of the equilibration of the polymer in the NpT ensemble. The longitudinal length of the simulation box in z-direction is calculated from the vapor-phase density of argon and the nominal thickness of the equilibrated polymer film h ) 10 nm. This gives for example Lz) 532.87 nm for an argon vapor density of 10 g/cm3 at TPE ) 80 K. We choose a large number of argon atoms in order to obtain a large box length perpendicular to the surface for a given density. A narrow vaporphase volume between the substrate surfaces would be similar to a vapor in a slit pore. In such confinement, the properties such as vapor pressure, density, and the critical properties differ from the corresponding bulk ones.21,22 Confinement can also affect phase transitions as condensation,21 freezing, and melting.23 A large vapor-phase box avoids these confinement effects.12,16 The states of the initial argon vapor are plotted in Figure 1 into the phase diagram of pure argon. It shows the saturated vapor densities (binodal) for argon obtained from experiments (Ar exp)24 and calculated from a Lennard-Jones equation of state (LJ-EOS).25 The phase diagram of the same potential model is required as used in the nucleation simulation in order to define the supersaturation of the vapor phase and also to establish a consistent comparison with the classical nucleation theory. For this purpose we calculated the phase equilibria and surface tension of the Lennard-Jones potential with a cutoff radius of 2.5 σ (LJ-2.5σ) by MD simulations. The corresponding data points are also added in Figure 1 together with a correlation curve. In order to locate the metastable region of the phase diagram the stability limit (spinodal) is required.

15786 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Rozas and Kraska

It can be estimated from a suitable equation of state applying the mechanical stability condition.26 Figure 1 shows that the state points used in the simulations here are all above the spinodal estimated from the LJ-EOS. We have no equation to estimate the spinodal for the LJ-2.5σ system, but from the comparison of different systems, it is known that the spinodal moves in the same direction as the binodal.26 It follows that the spinodal of the LJ-2.5σ system is very likely below the spinodal of the LJ-EOS and hence the simulation state points are clearly inside the metastable region of the pure argon phase diagram. This estimation does of course not include the substrate interaction. 2.2. Classical Nucleation Theory. In the case of a pure substance, vapor-phase nucleation is a process of the appearance of nanoscopically small molecular clusters of a new condensed phase by density fluctuations. The size of a cluster can be characterized by the number n of molecules constituting the cluster. The clusters of size n ) n* at the activation barrier of the nucleation process are called nuclei (or critical nuclei). The smaller (n < n*) or larger (n > n*) ones are called subnuclei or the supernuclei, respectively. The supernuclei can grow spontaneously to macroscopic sizes. Heterogeneous nucleation takes place when the old and forming phases are in contact with a third phase or molecular species. According to the shape and dimension of the clusters on a substrate different kinds of heterogeneous nucleation can be defined. Two idealizations for the vapor condensation on a substrate are perused here; clusters grow as droplets forming a contact angle with the surface (HEN 3D) and clusters grow as disks of fixed thickness along their periphery (HEN 2D). The work of formation of a three-dimensional condensed cluster of size n on a surface is given by the sum of the bulk contribution and the surface contribution:27

∆Gn ) - n∆µ + aγeffn2/3

(1)

The supersaturation S is the driving force for the nucleation. It is related to the chemical potential difference between the perfect gas vapor phase and the bulk liquid phase by ∆µ ) kBT ln(S). The surface term is the product of the cluster surface and the effective value of the specific surface energy. For cap-shaped clusters it is a ) (36πV02)1/3, with V0 the volume occupied by a molecule in the cluster, and further γeff ) ψ1/3(θ) γ, with ψ the Volmer function of the contact angle θ given by27

ψ(θ) ) (1/4)(2 + cos θ)(1 - cos θ)2, θ ∈ [0°,180°] (2) Here θ is the angle between the surface and the tangent plane on the cluster basis. When the relative interaction between a surface and the adatoms is strong enough in comparison to the interaction between adatoms, the contact angle defined by the Young relation loses its usual meaning and the two-dimensional nucleation model is more suitable. In contrast to the threedimensional growth here, only a part of the cluster surface is exposed to the vapor phase and can take part in the growth. Hence the monomers can attach only to a cluster surface area given by A ) bhn1/2, where b ) 2(πa0)1/2 is a shape factor for a disk-shaped cluster and h is the height of the cluster. The work of formation of a n-atom cluster is given by27

∆Gn ) - n(∆µ - a0∆γ) + bκn1/2

molecule in the cluster and ∆γ is a specific surface energy defined as

∆γ ) γ + γi - γs

(4)

Here γ, γi, and γs are the specific surface energy of the interfaces cluster-vapor, cluster-substrate, and substrate-vapor respectively. In terms of this parameter, we can define a crossover from the vapor homogeneous nucleation (HON), passing trough the growth of three-dimensional clusters on a substrate (HEN 3D), to the formation of layers on a substrate (HEN 2D). The conditions ∆γ > 0, ∆γ ) 0, and ∆γ < 0 correspond to incomplete, complete, and “better than complete” wetting respectively.27 The nucleation rate J, giving the number of supernuclei generated in the system per unit time and per unit volume (HON) or area (HEN), is an important kinetic characteristic of the nucleation process. A general classical expression for the stationary nucleation rate is given by the product of the Zeldovich’s z factor, the attachment frequency of monomers fn*, and the nuclei concentration Cn* in the stationary regime:27

J ) z fn*Cn*

(5)

The Zeldovich’s factor z and the stationary nuclei concentration Cn* depend on the expression of the formation work of a n-atom cluster ∆Gn for the different types of nucleation. The stationary nuclei concentration is given by the Boltzmann cluster size distribution27

Cn* ) C0 exp( - ∆G*/kBT)

(6)

The nucleation barrier energy, corresponding to the maximum of the function ∆Gn, for the three- and two-dimensional models are respectively, as follows:

∆G* )

3 3 b 2κ 2 4 a γeff and ∆G* ) 2 27 ∆µ ∆µ - a0∆γ

(7)

As shown by Zeldovich z is approximately given by z ) [-(d2∆Gn/dn2)n)n*/2πkBT]1/2 so that we obtain here

z ) ∆µ2/(64π2kBT)1/2V0γeff3/2 and z ) (∆µ - a0∆γ)/(4πkBT)1/2a01/2κ (8) for HEN 3D and for HEN 2D models respectively. The attachment frequency is a kinetic term that depends on the transport mechanism of monomers toward the growing cluster surface. Here we express the attachment frequency of the condensing vapor fn* by the product of the impingement flux I, the surface area A of the cluster where the monomers can attach, and the fraction of impinging monomers R, which effectively attach to the cluster surface. The impingement flux is approximated by the Hertz-Knudsen relation:

I)

p (2πmkBT)1/2

(9)

(3)

The island boundary contribution is defined in terms of a specific edge energy or line tension κ ≈ γh while the driving force includes an additional term -a0∆γ, where a0 is the area of a

Here p is the pressure of the vapor, m the mass of the atoms/ molecules, kB the Boltzmann’s constant, and T the temperature. This flux depends only on the temperature and pressure of the vapor phase. The following expression for the nucleation rate

Heterogenous Nucleation

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15787

is obtained in particular for cap-shaped clusters considering the contact area of the cap with the vapor:17

J ) RC0V0

( )( ) (

(1 - cos θ) 2γ 2 πmψ

1/2

)

p ∆G* exp kBT kBT

(10)

In the two-dimensional growth mechanism the monomers attach to the periphery of a cluster and the resulting expression for the nucleation rate derived here is

(

J ) RC0bh

)( ) (

∆µ - a0∆γ 8πm

1/2

)

p ∆G* exp kBT kBT

(11)

The concentration of nucleation sites C0 is approximated as the monomer density for the homogeneous nucleation C0 ) F. For HEN C0 is the concentration of active sites on the substrate surface here approximated as C0 ) F2/3. 2.3. Determination of Atomic Layers. The properties of argon in the vapor phase and on the polymer surface are calculated separately. The argon atoms adsorbed on the polyethylene surface are recognized with a three-step algorithm: (a) first the CH2 groups in the film surface are localized, (b) the argon atoms in the first adsorbed layer are determined, and (c) step b is repeated layer by layer until no more adsorbed atoms are found. The recognition of the surface CH2 groups is based on an implementation of the cone method.28 A polymer group site is located at the surface if a cone with its vertex located in the considered site and oriented perpendicular to the surface does not contain any other group. In the analysis, here the radius and angle of the cone are defined as 1 nm and 23° respectively. The first layer of adsorbed argon atoms is recognized by a search of all the argon atoms separated by a distance smaller than a certain cutoff from each surface CH2 group, which is the cluster definition by Stillinger’s criterion 1.5σ.29 In Figure 2c, the result of that algorithm for the first layer is shown. The search is repeated for the atoms deposited on higher layers as shown in Figure 2c. In order to recognize clusters on the surface the Stoddard algorithm30 is used. 3. Results In Figure 2a, a snapshot after 5 ns of a condensation simulation at TPE ) 80 K and F ) 50 g/dm3 is shown. A small region around the polymer film is exhibited only while the box perpendicular to the film is much larger, as mentioned above. At this temperature the polymer keeps to a large extent its crystalline structure as setup in the initial state. The surface of the polymer film (Figure 2b) is well ordered. At this point in time the first argon layer on the substrate is almost completely filled (Figure 2c). The higher layers are not filled, one rather observes surface clusters indicating deviations from layer-bylayer growth (Figure 2d). In order to illustrate the progression of the condensation process simulation with respect to the argon deposition and temporary temperature gradients in the system the density and temperature profiles of argon and the substrate are plotted in Figure 3. Two different points in time for a simulation with an initial vapor density of 50 g/dm3 and TPE ) 80 K are chosen. At 0.5 ns (Figure 3a), the density of the deposited film of argon atoms is lower than the saturated liquid density of LJ2.5σ argon at the same temperature. This is because the first argon layer is not completely filled yet. After 5 ns the argon film density reaches the value of saturated liquid LJ2.5σ argon at given temperature. In the first layer the density even slightly exceeds this value (Figure 3b) likely caused by surface induced close

Figure 2. Snapshots of a condensation simulation at TPE ) 80 K and F ) 50 g/dm3 after 5 ns. (a) Complete polymer film with condensed argon and a part of the argon vapor phase. (b) Surface layers of polyethylene. (c) Polyethylene surface layer and first argon layer. (d) Polyethylene surface layer with first and all subsequent argon layers. Each layer has a different grayscale (red shade).

packing. At 0.5 ns the system is not equilibrated yet and the vapor has developed a temperature gradient from the bulk vapor phase to the polymer surface at TPE ) 80 K (Figure 3a). This is because the thermostat is applied to the polymer only. The temperature at the polymer surface is slightly above that of the center of the film. After 5 ns, the temperatures of the vapor phase and of the film, including its surface layer, have approached (Figure 3b). There are large temperature fluctuations in the vapor-phase caused by the low density of the remaining vapor phase and hence the low number of atoms for averaging.

15788 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Figure 3. Density and temperature profiles at (a) 0.5 ns and (b) 5 ns for a simulation at TPE ) 80 K and F ) 50 g/dm.3

Another important property illustrating the actual course of the simulation is the change of the temperature during the condensation process as shown in Figure 4 for the different domains of the system. At low vapor-phase density (Figure 4a) in given time only a small amount of argon is deposited on the substrate. The condensed argon heats up the surface layer of the substrate by a few degrees only. The condensed argon film rapidly reaches a maximum temperature at the beginning of the condensation and then slowly decreases. Within about 5 ns, the temperature of the argon film reaches the value of the substrate. The temperature of the argon vapor phase starts at the initial

Rozas and Kraska temperature of 200 K and decreases very slowly due to the heat exchange with the substrate. At higher vapor-phase density (Figure 4b) we find similar temperatures for the substrate surface layer and the substrate bulk as for low vapor-phase density (Figure 4a). The temperature of the condensed argon also reaches the substrate temperature in about 5 ns. A major difference is observed for the vapor-phase temperature which starts here for all simulations at 200 K. The fast thermal equilibration of the vapor at higher density can be explained by the high number of collisions and evaporation of argon atoms from the condensed film. Also the effective expansion of the vapor-phase caused by the withdraw of atoms by condensation contribute to the temperature decrease: at the very beginning, the vapor behaves as a system expanded under quasi-adiabatic conditions because the decrease of internal energy caused by the condensation is faster than the energy exchange with the surface. This effect can be observed as a reduction of the vapor temperature below the temperature of the surface and becomes remarkable in simulations of vapors at high initial saturation (see Figure 4b at 7 ns) and also for a small number of vapor phase atoms. Avoiding such effect to a large extent is another reason why we here use the relatively large number of 5000 argon atoms. In order to analyze the growth of the first and all subsequent argon layers, we plot in Figure 5 the number of ad-atoms for each layer separately at different vapor-phase densities. At low vapor density (Figure 5a) in the beginning, e.g., the first 2 ns, most atoms are adsorbed in the first layer. The growth of the second and third layers follows slowly, while in given time almost no atoms are in the fourth and higher layers. This system is close to the binodal (see Figure 1), and since only the first layer is filled completely in the equilibrium state the deposition process is close to adsorption. With increasing density (Figure 5b) the growth becomes faster in the first layer and in the subsequent layers. Since more or less one layer fills after the other, we have here two-dimensional layer-by-layer growth. Scaling the time axis by the vapor-phase density does not give congruent plots as one might expect. For example, as a reference point, we may take the number of atoms in the third layer being approximately 90 at 12 ns and F ) 10 g/dm3 (Figure 5a). For a three times higher density, F ) 30 g/dm3, such scaling would suggest a three times faster growth and the reference point would correspond to 4 ns. However, we find approximately 580 in the first layer at 4 ns and F ) 30 g/dm3 and not 90. For a six times higher density, F ) 60 g/dm3, and one-sixth of the simulation time of the reference point (2 ns), we even find approximately 650 atoms in the third layer (Figure 5b).

Figure 4. Evolution of the temperature of the argon in the vapor phase, condensed argon, the surface layer of the system (atoms or sites in contact with vapor) and the polymer bulk film for (a) a high and (b) a low argon vapor density.

Heterogenous Nucleation

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15789

Figure 5. Argon ad-atoms per layer at TPE ) 80 K and two different argon vapor-phase densities.

Figure 6. Number of atoms distributed on the first and on the higher adsorbed layers as a function of the total number ad-atoms a for (a) TPE ) 80 K and (b) for TPE ) 60 K. Perfect layer-by-layer growth for the first layer is indicated by the dashed diagonal.

Therefore, a simple density scaling of the time axis is not possible here. In Figure 6, the number of ad-atoms in the different layers is plotted as a function of total number of adsorbed atoms for one simulation at TPE ) 80 K (Figure 6a) and one at TPE ) 60 K (Figure 6b). For the first layer perfect layer-by-layer growth is indicated by the dashed diagonal. One can recognize a deviation from layer-by-layer growth for the first layer with increasing density from 10 to 60 g/dm3 for 80 K and from 5 to 30 g/dm3 for 60 K, which are the density ranges investigated here limited by the width of the metastable region in the phase diagram (Figure 1). For comparison of parts a and b of Figure 6, one should keep in mind that the range of density and supersaturation is almost twice as large for 80 K. In the case of perfect layerby-layer growth the number of atoms in the subsequent layers (>first layer) should be zero until the first layer is completely filled, which corresponds to approximately 870 atoms at 80 K and 920 atoms at 60 K. The nominal areas of the substrate surface are slightly different for the simulations at 60 and 80 K. Hence the number of atoms per unit area is more suitable for comparison. It is 7.31 atoms/nm2 at 60 K and 6.87 atoms/ nm2 at 80 K. This difference of about 6.4% in the twodimensional concentration of the atoms in the first layer at different temperatures is related to the thermal expansion of the liquid. The bulk liquid densities Fl calculated here from simulations of equilibrated films for the LJ2.5σ argon are 1480.1 g/dm3 at 60 K and 1340.8 g/dm3 at 80 K. In order to compare the difference in the two-dimensional concentration of atoms on the substrate surface to that in the three-dimensional

bulk liquid density Fl, we have to scale the density as Fl2/3. This relation is exact for a perfect crystal and an approximation for the liquid here. It gives a difference of 5.8% in density, which is comparable to that of the concentration of atoms on the surface. From the simulation results the nucleation rates on the substrate surface are calculated by the method of Yasuoka and Matsumoto.31,32 This method is illustrated in an exemplary manner in Figure 7. It shows the evolution of the number of clusters larger than a given threshold value. In this diagram, plots for threshold values from 5 to 25 atoms are shown. The slope of the linear ascent gives the nucleation rate after division by the nominal surface area 2LxLy. Figure 7b shows that the slope converges to a limiting value with increasing threshold value. Such a decreasing slope corresponds to coalescence processes on the surface, which leads to a decrease of the number of clusters. Therefore, smaller threshold values correspond to an early stage of growth where nucleation dominates over growth by coalescence. Here we use the threshold value 5 (squares in Figure 8) for the calculation of the nucleation rate but also in the limit at 25 (circles in Figure 8) in order check the effect of coalescence on the results. The application of the method of Yasuoka and Matsumoto can be problematic in the extremes of low and high nucleation rates when it is applied to small finite surfaces. At high a nucleation rate, the regime of linear ascent in cluster growth exhibits a short duration due to the early coalescence of the clusters into a single layer. At low nucleation rates, the slope converges at a threshold size similar to or higher than the quantity of sites on the surface and therefore

15790 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Rozas and Kraska

Figure 7. (a) Number of cluster larger than a given threshold size as a function of the simulation time for an argon vapor density of 30 g/dm.3 The dashed lines are regression in the domain of linear ascend. (b) Dependence of the slopes of the linear domain plotted in part a on the threshold cluster size. The data point for 60 g/dm3 and n ) 25 is not available due to poor statistics.

Figure 8. (a) Comparison of the nucleation rates obtained from the simulations (symbols) with the heterogeneous three-dimensional classical nucleation model (HEN 3D) with two different contact angles and with the heterogeneous two-dimensional classical nucleation model (HEN 2D). Both models are applied at T ) 85 K. (b) Comparison of the nucleation rates with the heterogeneous two-dimensional classical nucleation model (HEN 2D). In both diagrams the circles are obtained from the limiting threshold value, the squares are obtained using the threshold value 5. The temperatures of the data points are 65 and 85 ( 2 K.

cannot be determined. Both difficulties can be avoided by increasing the size of the surface; however, that also increases the required computation power. The nucleation rates obtained from the simulation are plotted in Figure 8 and compared to different versions of the classical nucleation theory presented above. In order to establish a consistent comparison between theory and simulation for all equilibrium properties such as the liquid and vapor densities, vapor pressure, and the interfacial tension, the values obtained from equilibrium simulation with the LJ-2.5σ model are used. The supersaturation ratio is approximated by S ) F/Feq using bulk densities. The theory was locally applied in the region close to the surface where the vapor condenses. The temperature is calculated from the argon atoms and the polymer sites located in the surface of the film. In contrast to the classical Szilard model,33 the temperature of the surface as well as the density of the vapor phase change during the simulation. Within the Szilard model, the clusters grow according to a reversible bimolecular mechanism, the monomers attach to clusters of size j to build another of size j + 1 ignoring coalescence. Therefore, averaged values in the time interval of the linear ascent of the growth curves (Figure 7a) are used as input for the classical nucleation theory. As mentioned above the temperature of the argon vapor exceed that of the substrate. It is roughly 4-7 K above the substrate temperature. We chose the argon vapor

temperature because it determines the supersaturation and calculate the curves of the classical theory at 65 and 85 K, which roughly corresponds to the argon vapor-phase temperatures. Simulations results at TPE ) 80 K indicate that argon deposition continues for saturation ratios lower than 1, which is actually in the stable region of the phase diagram for pure argon. Here we have a surface present at which vapor atoms can adsorb even though the system is stable with respect to the pure fluid phase diagram. According to the observed growth mode and to the spatial distribution of the atoms on the surface, the HEN 3D model at high wettability (low contact angles) should represent the nucleation rates obtained from simulation. The results shown in the Figure 8 corroborates this observation. The range of application of the HEN 3D model becomes wider the lower the contact angle is. The HEN 2D and the HEN 3D (for low contact angles) models show similar behaviors at high supersaturation. It suggests that the nucleation rate changes are less sensitive to the growth mechanism at higher saturation. Nevertheless, at low saturation, the three-dimensional growth model predicts a fast decay of the nucleation rate while the simulation results maintain the trend of the high saturation regime. The HEN 2D model seems to be suitable to explain the nucleation or rather adsorption observed outside the metastable region of the pure argon phase diagram, S < 1, where the HEN 3D model

Heterogenous Nucleation

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15791

cannot be applied. The general distinctions between adsorption and heterogeneous nucleation are the saturation state of the vapor and the number of layers on the substrate. In the case of a stable monolayer or a layer with fewer atoms, it is rather an adsorption process. If more and more layers continuously are formed until the vapor is in a saturated state, it is a heterogeneous nucleation process. The specific surface energy ∆γ turns to be the relevant parameter to describe the trend of the J-S data. The dependence of the effective specific surface energy of clusters ∆γ on the temperature can be interpreted in terms of Dupre´ relation:27

∆γ ) 2γ - βa

(8)

Here βa denotes the specific adhesion energy of a cluster and the surface. The specific energy of the interface between liquid clusters and the vapor decreases at higher temperature while the adhesion energy remains almost constant. 4. Summary and Conclusion The investigation of heterogeneous nucleation for the weakly interacting polyethylene-argon system and with a relatively strong adhesion exhibits two-dimensional growth with an increasing trend to three-dimensional growth for rising supersaturation of the vapor. Even though the system exhibits good wettability, three-dimensional islands on the first absorbed layers are formed at higher supersaturation. We find that a simple scaling of the time axis by the density to obtain coinciding growth curves is not possible for heterogeneous nucleation as it is for homogeneous nucleation. This is in agreement with our earlier investigation on metal deposition on polyethylene.19 The use of a nonisothermal ensemble unveils temperature gradients in the system which are expected in an experiment. First, a temperature gradient is developed in the vapor because of the heath transfer between the vapor and the surface at lower temperature. Second, formation of a liquid argon film on the polymer substrate is accompanied by heat exchange between the vapor phase and the substrate. Both processes lead to an evolution of the temperature of the vapor, the substrate bulk and the surface, which have to be accounted for in the analysis of the simulation data, for example by classical nucleation theory. A comparison of the simulation results with classical nucleation theory shows that the heterogeneous nucleation rates can be modeled well by classical nucleation theory using the macroscopic properties of the same model fluid as used in the nonequilibrium simulations. The differences in the growth modes are less important in the regime of high saturation. Both versions of the classical nucleation theory: HEN 3D applied at low contact angle and the two-dimensional version of HEN 2D developed here, well describe the J-S data here obtained from the simulation. However, at low supersaturation the selection of the model is critical, due to the high wettability of this system the twodimensional heterogeneous version of the classical nucleation

theory is the most suitable model to describe the simulation data. While the three-dimensional model predicts the cessation of nucleation near the binodal curve of the pure vapor, S f 1, the two-dimensional growth model predicts a smooth transition from nucleation to adsorption of undersaturated vapor. Acknowledgment. R.R. has been supported by a Ph.D. fellowship of the Deutscher Akademischer Austauschdienst (DAAD). This work has been partially supported by DFG Project KR1598/19-2. Furthermore, the Regionales Rechenzentrum at the University of Cologne is acknowledged for providing computing resources and support. References and Notes (1) Ravishankar, N.; Carter, C. B. Philos. Mag. Lett. 2005, 85, 523531. (2) DeMott, P. J.; Cziczo, D. J.; Prenni, A. J.; Murphy, D. M.; Kreidenweis, S. M.; Thomson, D. S.; Borys, R.; Rogers, D. C. Proc. Nat. Acad. Sci. U.S.A. 2003, 100, 14655-14660. (3) Becker, R.; Do¨ring, W. Ann. Phys. 1935, 24, 719-752. (4) Volmer, M.; Weber, A. Z. Z. Phys. Chem. (Leipzig) 1926, 119, 277-282. (5) Volmer, M. Z. Electrochem. 1929, 35, 555-561. (6) Zeng, X. C.; Oxtoby, D. W. J. Chem. Phys. 1991, 94, 4472-4478. (7) Padilla, K.; Talanquer, V. J. Chem. Phys. 2001, 114, 1319-1325. (8) Yasuoka, K.; Matsumoto, M. J. Chem. Phys. 1998, 109, 84518462. (9) Laasonen, K.; Wonczak, S.; Strey, R.; Laaksonen, A. J. Chem. Phys. 2000, 113, 9741-9747. (10) Toxvaerd, S. J. Phys. Chem. 2001, 115, 8913-8920. (11) Lu¨mmen, N.; Kraska, T. J. Aerosol Sci. 2005, 36, 1409-1426. (12) Kholmurodov, K. T.; Yasuoka, K.; Zeng, X. C. J. Chem. Phys. 2001, 114, 9578-9584. (13) Senger, B.; Schaaf, P.; Corti, D. S.; Bowles, R.; Pointu, D.; Voegel, J.-C.; Reiss, H. J. Chem. Phys. 1999, 110, 6438-6450. (14) Chen, B.; Siepmann, J. I.; Oh, K. J.; Klein, M. L. J. Chem. Phys. 2002, 116, 4317-4329. (15) Zapadinsky, E.; Lauri, A.; Kulmala, M. J. Chem. Phys. 2005, 122, 114709. (16) Toxvaerd, S. Mol. Simul. 2004, 30, 179-182. (17) Kimura, T.; Maruyama, S. Microscale Thermophys. Eng. 2002, 6, 3-13. (18) Maruyama, S. Microscale Thermophys. Eng. 1998, 2, 49-62. (19) Rozas, R.; Kraska, T. Nanotechnology 2007, 18, 165706. (20) Sumpter, B. G.; Noid, D. W.; Wunderlich, B. Macromolecules 1992, 25, 7247-7255. (21) Thommes, M.; Findenegg, G. H. Langmuir 1994, 10, 4270-4277. (22) Jiang, S.; Rhykerd, C. L.; Gubbins, K. E. Mol. Phys. 1993, 79, 373-391. (23) Alba-Simioesco, C.; Coasne, B.; Dosseh, G.; Dudziak, G.; Gubbins, K. E.; Radhakrishnan, Sliminska-Bartkowiak, M. J. Phys.: Condens. Matter 2006, 18, R15-R68. (24) IUPAC. International Thermodynamic Tables of the Fluid State Argon; Butterworths: London, 1972. (25) Kolafa, J.; Nezbeda, I. Fluid Phase Equilib. 1994, 100, 1-34. (26) Kraska, T. Ind. Eng. Chem. Res. 2004, 43, 237-242. (27) Kashchiev, D. Nucleation; Butterworth-Heinemann: London, 2000. (28) Wang, Y.; Teitel, S.; Dellago, C. J. Chem. Phys. 2005, 122, 214722. (29) Stillinger, F. H. J. Chem. Phys. 1963, 38, 1486-1494. (30) Stoddard, S. D. J. Comput. Phys. 1978, 27, 291-293. (31) Yasuoka, K.; Matsumoto, M. J. Chem. Phys. 1998, 109, 84518462. (32) Yasuoka, K.; Matsumoto, M. J. Chem. Phys. 1998, 109, 84638470. (33) Schmelzer, J.; Ro¨pke, G.; Mahnke, R. Aggregation Phenomena in Complex Systems; Wiley-VCH Verlag: Weinheim, Germany, 1999.