2360
J. Phys. Chem. C 2009, 113, 2360–2367
Accommodation of Gases at Rough Surfaces Natasa Mateljevic, Jay Kerwin, Sharani Roy, J. R. Schmidt, and John C. Tully* Department of Chemistry, Yale UniVersity, New HaVen, Connecticut 06520-8107 ReceiVed: September 1, 2008; ReVised Manuscript ReceiVed: December 3, 2008
A simple numerical model is proposed to compute the energy and momentum accommodation of molecules scattered from highly corrugated, disordered surfaces. The model is an extension of the “washboard model”, which assumes that the component of the molecule’s momentum parallel to the local surface tangent is conserved on impact and the normal component is altered by a hard, elastic collision with a moving surface “cube” with an adjustable effective mass. The surface is represented by Gaussian hills and valleys of random location and height. In contrast to the washboard model, the current model is fully three-dimensional and includes in-plane and out-of-plane scattering as well as trapping-desorption. In addition, it can be applied to highly corrugated surfaces and does not invoke a regular, periodic topography. This increased realism comes at the expense of an analytical solution; numerical simulations must be performed. We develop a very efficient procedure for carrying out the simulations. We test the model by comparing detailed angular and velocity scattering distributions for Xe scattering from a Pt(111) surface with those obtained by realistic molecular dynamics simulations of the same system. We then apply the model to the accommodation of H2, N2, CO, and CO2 on surface materials employed on vehicles in low Earth orbit. The model is capable of accurately reproducing the results of experimental measurements on these highly corrugated surfaces. I. Introduction The transfer of energy and momentum resulting from the collision of a gas molecule with a surface has been a subject of fundamental interest for a many years.1 Moreover, knowledge of energy and momentum accommodation is required to understand practical phenomena, such as the drag and frictional heating of an object moving through a rarefied gas, adsorption and chemical reaction at surfaces, erosion of surfaces of reentry vehicles and satellites in low Earth orbit, and hydrodynamic modeling of gases in contact with surfaces. Our understanding of gas-surface energy transfer has evolved from the rudimentary “hard cube model”1-4 to our current ability to carry out three-dimensional molecular dynamics simulations with accurate force fields,5 with hundreds of surface atoms to properly represent phonon motions, and even with inclusion of electronic dissipation at metal surfaces.6 Laboratory studies have been revolutionized by molecular beam technology and the myriad of molecular and surface characterization tools now available. Both molecular beam and molecular dynamics studies are demanding, however, and the latter require detailed prior knowledge of the system. There remains a need for relatively simple models that can provide scattering distributions and accommodation coefficients of useful accuracy for relatively uncharacterized systems. This is especially true for highly corrugated or disordered surfaces that are commonly used in practice and for which current models are generally inadequate. We present here a simple numerical model applicable to rough, disordered surfaces. The model is an outgrowth of the “washboard model”,7 which in turn is an extension of the “hard cube model”.2,3 The hard cube model was proposed to describe the scattering of an atom from a very smooth surface. The basic idea of the hard cube model is simple: if the surface is smooth, the colliding atom will experience only very small forces in the plane of the surface. Thus, the in-plane momentum is * Corresponding author. E-mail:
[email protected].
assumed to be unchanged. The component of gas momentum normal to the surface is altered through a hard collision between objects of mass m (gas) and M (effective surface “cube”), with the initial velocity of the surface cube chosen at random from a Maxwellian distribution at the surface temperature T. The only parameter required by the model, in addition to the experimentally defined properties T, m, the incident gas energy, Ei, and the angle of approach, θi, is the effective mass, M, of the surface cube. This parameter regulates the overall extent of inelasticity; that is, the “softness” of the surface. A modification of the hard cube model (sometimes called the “soft cube model”)1,4 introduces an additional parameter, an attractive potential well of depth W, to incorporate trapping (sticking). The hard and soft cube models have proved surprisingly successful in qualitatively reproducing observed trends in nonreactive scattering of gases from smooth, single-crystal surfaces. However, they are limited by the assumption that the gas-surface interaction is ideally flat; the atom experiences no in-plane forces. A number of models have been subsequently proposed that incorporate surface corrugation.7-13 One of these, the “washboard model”7, is a direct extension of the hard cube model. The washboard model retains the hard cube assumptions that the lateral component of momentum is conserved and the normal component is modified via a hard (elastic) collision. However, in the washboard model, the normal and lateral directions are defined locally by the slope of the surface at the point of impact of the molecule. The surface is assumed to exhibit a sinesoidal profile with a maximum slope A. Thus, the washboard model invokes three adjustable parameters, M, W, and A, in addition to the experimentally assigned properties. In cases in which the corrugation parameter is sufficiently small that multiple scattering events and “shadowing” are not significant, an analytical expression has been derived for the washboard model angular and velocity scattering distributions, energy accommodation coefficient, and trapping probability.7 The model has been shown to be a significant improvement
10.1021/jp8077634 CCC: $40.75 2009 American Chemical Society Published on Web 01/20/2009
Accommodation of Gases at Rough Surfaces
J. Phys. Chem. C, Vol. 113, No. 6, 2009 2361
over the hard cube model, even for scattering from “smooth” closely packed faces of metallic crystals.13-16 However, it is not adequate for highly corrugated surfaces that exhibit shadowing, multiple scattering, and significant out-of-plane scattering. Recently a “washboard with moment of inertia” (WMI) model was proposed to overcome some of the deficiencies of the original washboard model.17 The WMI bestows on the effective surface object a moment of inertia in addition to its effective mass, M. This introduces an anisotropic response of the surface, a crucial effect observed in molecular dynamics simulations of atomic scattering from highly corrugated organic polymer surfaces.18 The long-chain polymer molecules are stiffer when impacted head-on and floppier when impacted from the side; the chains can wag laterally relatively freely compared to compression of the chain length. The WMI model was quite successful in reproducing the 2D energy-angle distributions for neon scattered from a decane thiol monolayer adsorbed on Au(111), as computed by molecular dynamics simulation.17,18 The WMI model, however, like the washboard model, assumes a periodic two-dimensional sinesoidal surface topography; the model is not applicable to disordered surfaces, and it does not describe out-of-plane scattering, a dominant feature in scattering from highly corrugated surfaces. In this paper, we present an alternative variation of the washboard model that is fully threedimensional with a random surface topography. The model is based on classical mechanics and invokes the washboard model assumption of conservation of local tangential momentum, but applied to a 3-dimensional, random surface topography. Unlike the hard cube and washboard models, the current model does not appear amenable to closed-form solution. Accommodation coefficients, trapping probablilities, and scattering distributions must be obtained by numerical integration. We present an algorithm for carrying out the required numerical integrations very efficiently, requiring less than a microsecond per trajectory on a 3.5 GHz pentium processor, thereby allowing rapid and thorough exploration of trends and parameter space. The model is formulated in Section II, with some details relegated to Appendix A. In Section III, we test the assumptions of the model by comparison with angular and velocity distributions for xenon scattered from the (111) face of platinum, as computed by realistic molecular dynamics simulations. In Section IV, we apply the model to analyze a systematic body of experimental results on momentum and energy accommodation of the molecules H2, N2, CO, and CO2 scattered from four different “practical” surfaces of relevance to the space shuttle. The model is very successful in reproducing the experimental observations for this set of highly corrugated and disordered surfaces. Conclusions are presented in Section V, and details of the molecular dynamics simulations are given in Appendix B. II. The Model A. Kinematics. Let m, Ei, and θi be the gas molecule mass, incident kinetic energy, and incident angle with respect to the surface normal, respectively. The initial components of velocity of the molecule in the laboratory frame are:
Vy ) 0 Vz ) - (2Ei ⁄ m) cosθi
of generality in restricting the initial velocity vector to the x-z plane, since the surface is assumed to be random with no preferred directionality in the x and y directions. The z component of velocity increases in absolute value to V˜ z as the molecule enters a square well of depth W.
V˜ z ) - (V2z + 2W ⁄ m)1⁄2
(2)
Let R and β be the polar and azimuthal angles of the vector normal to the local tangent of the surface at the position that the trajectory intersects the surface. Only the polar angle R is illustrated in the two-dimensional plot of Figure 1. We transform the initial velocity to a frame in which the z axis is taken along the local surface normal direction:
ux ) cos R cos βVx - sin RV˜ z uy ) -sin βVx uz ) sin R cos βVx + cos RV˜ z
(3)
The fundamental approximation of the model is a threedimensional generalization of the hard cube and washboard model assumptions: in-plane momenta are conserved on impact, and the change in the normal momentum is governed by a onedimensional impulsive collision between objects with mass ratio µ ) m/M, where M is the effective mass of the surface cube.
ux ) ux uy ) uy uz )
( µµ +- 11 )u + ( µ +2 1 )V z
c
(4)
In eq 4, Vc is the initial thermal velocity of the surface object, selected from a Maxwellian distribution at the surface temperature, T. Transforming back to laboratory coordinates, the components of velocity after collision are Vx ) (1 - sin2 R cos2 β)Vx - cos R sin R cos βV˜ z + sin R cos βuz Vy ) -sin2 R sin β cos βVx - cos R sin R sin βV˜ z + sin R sin βuz
(5)
V˜ z ) -cos R sin R cos βVx + sin2 RV˜ z + cos Ruz
Accounting for the deceleration in escaping the well, the final z-component of velocity is
Vx ) (2Ei ⁄ m)1⁄2sinθi 1⁄2
Figure 1. Schematic illustration of a 2D off-normal incident trajectory and the necessary quantities to reweight PA(θi, R, β). The global surface normal is given by the vector Nˆ, and the local surface normal by the ˆ . The path of the incident trajectory is given by the vector Iˆ. vector Q
(1)
where z is the direction normal to the plane of the surface, which we will refer to as the global normal. Note that there is no loss
1⁄2 Vz ) (V˜ 2 z - 2W ⁄ m)
(6)
For the incident molecule to strike the surface cube and recoil with final z-component of velocity, V′z, sufficient to escape the attractive well, the following condition must be satisfied:
2362 J. Phys. Chem. C, Vol. 113, No. 6, 2009
V˜ z > (2W ⁄ m)1⁄2
Mateljevic et al.
(7)
If the condition of eq 7 is not satisfied, then the molecule is designated as trapped. Equations 1–7 define the kinematics of an individual trajectory that strikes the surface at a point where the direction of the local surface normal is defined by the angles R and β. To complete the model, a surface profile must be constructed with the desired topography and corrugation amplitude. A Monte Carlo sampling of trajectories could then be carried out, choosing the initial x and y positions of the gas molecule at random from uniform distributions and the velocity, Vc, of the surface cube at random from a Maxwellian distribution. This is feasible but tedious, primarily because the possibility of shadowing requires that each striking point be determined by a numerical search. It is much more efficient to first construct the probability distribution function PA(θi; R, β) that a molecule with incident angle θi will strike a surface of corrugation amplitude A (to be defined below) at a location where the polar and azimuthal angles defining the local surface normal are R and β. Once this distribution has been computed, all final scattering properties can be obtained simply by cycling through the relevant ranges of R, β, and cube velocity Vc, at sufficiently small increments, weighting each trajectory by the product PA(θi; R, β)PT(Vc)S, and summing over all trajectories, normalized by the total weight of all trajectories sampled. PT(Vc) is the unnormalized Maxwell-Boltzmann probability distribution function that at temperature T the surface cube will have velocity Vc along the local normal direction, and S is the relative collision velocity.
PT(Vc) ) exp(-MVc2 ⁄ 2kT)
(8)
S ) Vc - uz
(9)
B. Probability Distribution Function PA(θi; r, β). Any probability distribution PA(θi; R, β) that adequately represents the profile of an actual surface of interest can be employed. Consistent with the simplicity of the model and the usual paucity of detailed information about the atomic-level corrugation of the gas-surface interaction potential for most systems, we choose PA(θi; R, β) from a Gaussian random surface profile. Specifically, the surface is constructed from a sum of randomly positioned Gaussian hills and craters. Each Gaussian i has the same width, σ, a height (or depth) Ci chosen at random from a uniform distribution, -C < Ci < C, and location (xi, yi) chosen uniformly in the ranges -L < xi < L and -L < yi < L. As the number, N, of Gaussians per unit area increases, the resulting corrugation increases without bound. However, after normalizing the height at each point by N-1/2, the statistical properties of the surface become invariant for large N. We derive in Appendix A a simple expression for the average value of the tangent of the polar angle R in the limit of large N.
A ) 〈tan R〉 )
C π σ 24
1⁄2
( )
(10)
The corrugation parameter, A, controls the smoothness or roughness of the surface and can be chosen to achieve any desired corrugation amplitude. The corrugation parameter, A; the effective mass, M, of the surface object; and the depth, W, of the (optional) attractive well are the only adjustable parameters of the model. We show in Appendix A that, for a given choice of corrugation parameter A and initial incident angle θi, the probability distribution, PA(θi; R, β), can be accurately approximated by
PA(θi ;R, β) )
[
]
sin R π exp - 2 (sec2 R - 1) 2 3 4A cos R 4A × (sin θi tan R cos β + cos θi)
(11)
where PA(θi; R, β) is set equal to zero if the expression in eq 11 is negative, as occurs for angles that are geometrically inaccessible. All of the calculations reported in this paper employ the distribution function of eq 11. C. Calculation of Observables. We consider directed impact (fast projectile or molecular beam) conditions for which the initial molecular kinetic energy, Ei, and incident angle, θi, are fixed. Results for a thermal distribution of incident angles and kinetic energies can be obtained by a Boltzmann-weighted sum of results from specified incident angles and energies, or else directly by Monte Carlo sampling of initial angles and energies if information about specific initial conditions is not required. For the selected values of Ei and θi, as well as the parameters A, M, and W, we numerically integrate over R, β, and cube velocity Vc. If the trajectory satisfies the condition of eq 7, we compute the final x, y, and z velocities: V′x, V′y, and V′z, from eqs 1–6. From these, we calculate the final polar and azimuthal scattering angles, θF and φF; energy, EF; and speed, VF, from 2 1⁄2 θF ) tan-1((V2 x + Vy ) ⁄ Vz)
(12)
φF ) tan-1(Vy ⁄ Vx)
(13)
2 2 1⁄2 V2 x + Vy + Vz
VF ) (
)
1 EF ) mVF2 2
(14) (15)
We then construct any desired histograms or averages of these quantities, with each trajectory weighted by PA(θi; R, β)PT(Vc)S/ PTOT, where PTOT is the normalization denominator; that is, the sum of PA(θi; R, β)PT(Vc)S over all R, β, and Vc. To account for trapping eventsstrajectories that do not satisfy eq 7swe assume that all trapped trajectories completely accommodate to the surface temperature. We then add to the histograms the properly weighted contributions from trapping, assuming a cos θ angular distribution and a velocity distribution given by
PBoltz(V) dV ) V3 exp(-mV2 ⁄ 2kT) dV
(16)
The energy accommodation coefficient, ε, is defined as
ε)
Ei - < EF > Ei - 2kT
(17)
Note that if all trajectories trapped, the mean final energy would equal 2kT and the energy accommodation coefficient would be unity (full accommodation). However, the energy accommodation coefficient is not simply equal to the trapping probability because directly scattered molecules have altered kinetic energies. Normal and tangential momentum accommodation coefficients can also be extracted:
|Vz|- ηn ) |Vz|- (kT ⁄ m)1⁄2 ηt )
Vx - Vx
(18) (19)
Finally, following Cook et al.,19,20 we define normal and tangential “reduced force coefficients” fn and ft:
Accommodation of Gases at Rough Surfaces
fn ) ft )
|Vjz| + |Vjz|
(Vx2 + Vy2 + Vz2)1⁄2 |Vjx| - |Vjx|
(Vx2 + Vy2 + Vz2)1⁄2
J. Phys. Chem. C, Vol. 113, No. 6, 2009 2363
(20)
(21)
Equations 20 and 21 avoid the singularities of eqs 18 and 19. III. Comparison with Molecular Dynamics Simulations We carried out a series of molecular dynamics (MD) simulations of the scattering of xenon atoms from the (111) face of a platinum crystal. Details of the MD simulations are given in Appendix B. The interaction potentials employed are approximate but are sufficiently realistic to provide a good test of the model for a case in which all interaction parameters are known. It is admittedly curious to apply the random corrugation model that is intended for irregular and highly corrugated surfaces to the very smooth and regular Pt(111) surface. Nevertheless, this test demonstrates that the model is useful even in the small corrugation limit. More importantly, because all MD interactions have been specified, this comparison serves to evaluate whether the optimized parameters of the random corrugation model actually reflect the properties of the gassurface scattering system. The random corrugation model requires specification of three adjustable parameters: the mean corrugation amplitude, A; the effective surface mass, M; and an optional attractive well depth, W. The effective surface mass M governs the inelasticity of the surface; the average amount of energy deposited into the surface will be higher for small M and lower for large M. An attractive square well of depth W can be included to describe more realistically the trapping of molecules on the surface. The front edge of the square well is assumed to be parallel to the surface plane and beyond the highest feature of the surface. Thus, the spatial extent of the well is irrelevant, and only the well depth is needed. The well depth can be set equal to zero to eliminate one parameter, if desired. In addition to these adjustable parameters, the model requires specification of the following experimental variables: the surface temperature, T; the incident kinetic energy of the molecule, Ei; its incident angle with respect to the surface normal direction, θi; and the mass of the molecule, m. Figure 2 shows the final in-plane angular and velocity (speed) distributions of Xe scattered from Pt(111) computed by MD and by the random corrugation model are in excellent agreement. The parameters of the random corrugation model were chosen to be the same for all the results in Figure 2: A ) 0.06, M ) 600 amu, and W ) 13.35 kJ/mol. A ) 0.06 corresponds to an average surface normal, 〈R〉 ) 3.43° and reflects the very small corrugation of the smooth Pt(111) surface. M ) 600 amu is a reasonable value for the effective mass, about 3 times the mass of an individual platinum atom. An incoming xenon atom typically interacts with more than one surface atom, and the surface atom response is stiffened due to harmonic forces between platinum atoms. An effective mass of 2 or 3 times the mass of an individual surface atom is typically found in applications of the hard cube and washboard models. The attractive well depth of 13.35 kJ/mol is one-half of the isosteric heat of adsorption of xenon on Pt(111) at zero coverage, 26.7 kJ/mol. This is also quite reasonable: the incoming atom does not usually visit the minimum energy site of the interaction potential. The results were relatively insensitive to the value of W, as excellent fits were obtained for values from 40 to 100%
Figure 2. Comparison between molecular dynamics simulations (points) and the random corrugation model (curves) for scattering of xenon atoms from the Pt(111) surface. The incident angle of approach with respect to the surface normal is 45°, and the surface temperature is 300 K for all cases. The incident Xe kinetic energy is listed in each panel. The left-hand panels show the in-plane angular scattering distributions. The right-hand panels show the in-plane velocity distributions. The random corrugation model parameters were chosen to be the same for each case: A ) 0.06, M ) 600 amu, W ) 13.35 kJ/mol.
of the isosteric heat of adsorption. Below 40%, the bimodal structure of the fits started to disappear, and only direct scattering contributed to the final distributions. IV. H2, N2, CO, and CO2 Scattered from Solar Panel Array Material Cook et al.19,20 have carried out a systematic experimental study of momentum accommodation in the scattering of the molecules H2, N2, CO, and CO2 from a series of “real” surfaces that have been used on solar panel arrays in low Earth orbit. Using a torsion balance apparatus, Cook et al. measured both the magnitude and direction of the forces exerted on the surface by molecules from a supersonic nozzle incident at a chosen angle and velocity. The surfaces studied were a polyimide plastic with the trade name Kapton; SiO2-coated Kapton; solar panel material that was essentially SiO2; and aluminum coated with a thermal insulating material called Z-93. The gases H2, N2, CO, and CO2 are produced by the control jets on the Space Shuttle via the reaction of monomethyl hydrazine with nitrogen tetroxide. The experiments of Cook et al. were directed at characterizing the interactions of these molecules with nearby surfaces. This body of experimental results provides an excellent opportunity to apply the random corrugation model to scattering from surfaces that are highly corrugated, irregular, and not wellcharacterized. Cook et al. report their data as reduced force coefficients, eqs 20 and 21. Figures 3-6 show comparisons between the experimental results and those of the random corrugation model. For all these comparisons, we chose the attractive well depth, W, to be zero to reduce the number of fitting parameters. For each type of surface, we optimized the two remaining parameters, A and M. Note that we used the same values of A and M for every gas scattering from a given surface. Although corrugation is a property of both the projectile and the surface, we chose to use the same parameters for each gas
2364 J. Phys. Chem. C, Vol. 113, No. 6, 2009
Mateljevic et al.
Figure 3. Reduced force coefficients for H2, N2, CO, and CO2 incident on Kapton. Surface temperature is 300 K. Points are experimental data from Cook et al.19 Curves are results of random corrugation model, with corrugation parameter A ) 0.92 and effective mass M ) 44 amu. The well depth is taken to be zero. The normal and tangential components of the reduced force coefficients are denoted by squares and triangles, respectively.
Figure 5. Reduced force coefficients for H2, N2, CO, and CO2 incident on solar panel (SiO2) material. Surface temperature is 300 K. Points are experimental data from Cook et al.19 Curves are results of random corrugation model, with corrugation parameter A ) 0.72 and effective mass M ) 46 amu. The well depth is taken to be zero. The normal and tangential components of the reduced force coefficients are denoted by squares and triangles, respectively.
Figure 4. Reduced force coefficients for H2, N2, CO, and CO2 incident on SiO2-coated Kapton. Surface temperature is 300 K. Points are experimental data from Cook et al.19 Curves are results of random corrugation model, with corrugation parameter A ) 0.7 and effective mass M ) 46 amu. The well depth is taken to be zero. The normal and tangential components of the reduced force coefficients are denoted by squares and triangles, respectively.
Figure 6. Reduced force coefficients for H2, N2, CO, and CO2 incident on Z-93-coated aluminum. Surface temperature is 300 K. Points are experimental data from Cook et al.19 Curves are results of random corrugation model, with corrugation parameter A ) 1.45 and effective mass M ) 8 amu. The well depth is taken to be zero. The normal and tangential components of the reduced force coefficients are denoted by squares and triangles, respectively.
because N2, CO, and CO2 have similar molecular sizes, and all four molecules have very weak attractive interactions with these surfaces. While adjusting parameters separately for each gas would have produced somewhat improved fits, Figures 3–6 show quite acceptable agreement. The overall shapes of the reduced force coefficients are not only in accord with experiment, but also in quite good quantitative agreement, including the values of the maxima and the positions where the normal and tangential curves cross. Cook et al.19 report scanning electron microscope images of the four surfaces employed in their scattering experiments. Kapton, SiO2-coated Kapton, and solar panel material were all
smooth on a similar scale; Z-93-coated aluminum was rough on a scale 10 times longer.19 With this in mind, it might be expected that Kapton, SiO2-coated Kapton, and solar panel material have similar reduced force coefficients. This expectation is borne out by both the experimental and simulation data. Figures 3, 4, and 5 are very similar for H2, N2, CO, and CO2. The extracted simulation parameters A and M are also nearly the same, with the best fit parameters for Kapton showing a slightly higher effective corrugation. The results for Z-93-coated aluminum are unique. When an incident beam undergoes complete thermal accommodation at the surface, the incident molecules trap, lose memory of their
Accommodation of Gases at Rough Surfaces
J. Phys. Chem. C, Vol. 113, No. 6, 2009 2365
Figure 7. Polar plots of the angular distributions of scattering as predicted by the model. The polar angle θF is given by the radial distance from the origin and ranges from 0 to 90°. The azimuthal angle, φF, varies from 0 to 360° counterclockwise around the circle. The particles are incident from the left at 30° from the normal. Thus, an azimuthal angle of 0 corresponds to forward in-plane scattering, and an azimuthal angle of 180° is backward in-plane-scattering. (a) Xe scattered from Pt(111) with incident energy of 20 kJ/mol and all parameters taken to be the same as in Figure 2. (b) N2 scattered from solar panel material, with incident energy of 49 kJ/mol and all parameters the same as for Figure 5. (c) N2 scattered from Z-93, with incident energy of 49 kJ/mol and all parameters the same as for Figure 6.
initial conditions, and desorb from the surface equilibrated at the surface temperature. This is termed trapping-desorption. In eqs 20 and 21, the final normal and tangential components of velocity become the thermal average quantities: (πkT/2m)1/2 and zero, respectively.20 The incident velocity normal and tangential components can be represented by Vi cos θi and Vi sin θi, respectively, where(Vx2 + Vy2 + Vz2)1/2. This makes the reduced force coefficients, in the limit of total thermal accommodation
f thermal ) cos θi + √πkT ⁄ 2mVi2 n
(22)
f tthermal ) sin θi
(23)
In Figure 6, both the experimental and simulation data follow eqs 22 and 23 very closely for the plots of N2, CO, and CO2. This is due to nearly complete thermal accommodation of these incident species on the rough Z-93 surface. The experiment and simulation depart from eqs 22 and 23 in the case of H2 because H2 is lighter than the recoil mass, thus reducing the trappingdesorption and pushing the scattering toward the impulsive regime. V. Summary We have proposed a simple numerical model to compute the energy and momentum accommodation of molecules scattered from highly corrugated, disordered surfaces. The model is an extension of the “washboard model” which assumes that the component of momentum parallel to the local surface tangent is conserved on impact, and the normal component is altered via a hard-sphere collision. The surface is represented by Gaussian hills and valleys of random location and height. In contrast to the washboard model, the current model is fully threedimensional and includes in-plane and out-of-plane scattering as well as trapping-desorption. In addition, it can be applied to highly corrugated surfaces and does not invoke a regular, periodic topography. This increased realism comes at the expense of an analytical solution; numerical simulations must be performed. We have derived an accurate analytical expression, eq 11, for the probability of striking the surface at a position where the polar and azimuthal angles between the local normal and the surface normal are given by R and β, respectively. The final scattering distributions are obtained efficiently by numerically integrating over R and β as well as over the one-dimensional Maxwellian velocity distribution of
the surface object. As demonstrated above, the model successfully reproduces detailed in-plane angular and velocity scattering distributions for Xe scattering from a Pt(111) surface computed by realistic molecular dynamics simulations. In addition, the model accurately reproduces experimentally determined momentum accommodation coefficients of H2, N2, CO, and CO2 on surface materials employed on vehicles in low Earth orbit. The model requires specification of three adjustable parameters: M, W, and A. The well depth, W, can often be obtained from measured binding energies or estimated from calculated or assumed interaction potentials. If the surface is ordered, the corrugation parameter, A, may also be obtained from electronic structure calculations or known atomic radii. For disordered surfaces, topographies obtained by scanning probe microscopy may be useful for estimating A. It should be emphasized, however, that the corrugation is a property of both the surface and the probe; if the scanning probe tip is much larger than a molecular size, the measured corrugation will be correspondingly reduced. Nevertheless, we expect that there will be a strong correlation between the corrugation parameter, A, and the topography determined by scanning probe microscopy. Finally, the effective mass parameter, M, or more properly the mass ratio m/M, is a dynamical property that accounts for the response of the surface to impact by the molecule; that is, the “softness” or inelasticity of the surface. This cannot be obtained from static measurements or electronic structure calculations. For metal surfaces, M is typically found to be 2 or 3 times the mass of a single surface atom. For more complicated materials, such as the satellite surface materials studied here, further experience is needed before M can be estimated in advance with any confidence. An important feature of the model is that it is fully 3-dimensional; that is, it predicts out-of-plane scattering, expected to be very significant for highly corrugated surfaces. Examples of calculated (θF, φF) angular distributions are shown in Figure 7. Figure 7a shows that, for the very smooth Xe-Pt(111) interaction, the scattering intensity peaks near the φF ) 0 forward, in-plane direction, as expected. However, for scattering of N2 from the solar panel material, Figure 7b, the spread in the azimuthal angular distribution actually exceeds that of the polar angle θF; out-of-plane scattering is dominant. Finally, for N2 scattering from Z-93, for which the calculated trapping probability is nearly unity, Figure 7c shows that the angular distribution is isotropic in φF (essentially cos θF).
2366 J. Phys. Chem. C, Vol. 113, No. 6, 2009
Mateljevic et al.
The model presented here is simple enough to be easily and widely used to characterize the corrugation or “roughness” of surfaces from molecular beam scattering data, not only from crystalline and disordered solid surfaces, but also from liquid surfaces for which a large amount of data has appeared recently.21 In addition, the model can be incorporated into largescale computer codes that simulate friction and drag via hydrodynamic equations. Such codes currently use oversimplified assumptions about gas-surface interactions, such as “slip” and “stick” boundary conditions. A realistic treatment of energy and momentum accommodation that accounts correctly for surface roughness should greatly improve the accuracy of these simulations.
〈S′x2 〉 ) 〈S′x2 〉 )
〈
N
∑
-xi -yi 1 -xiCie 2σ2 6 2πFσ i)1 2
2
N
∑
-xjCje
-xj2-yj2 2σ2
j)1
〈
-xi N 〈Ci2〉 xi2e σ2 6 2πFσ
2
〉
)
C 〉〈e 〉 ) 12σ -yi2
2
σ2
2
(A5)
where we have again used the property that the random variables are uncorrelated. By the central limit theorem, we see that the surface is isotropic and its gradient is normally distributed with zero mean and variance C2/12σ2. Let the global normal Nˆ ≡ 〈0, 0, 1〉, the local surface normal b ≡ 〈Sx′, Sy′, 1〉, and the angle between at some position x, y Q the two be R:
b ) |Nˆ||Q b |cos R NˆQ Appendix A. The Probability Distribution PA(θi; r, β) Consider a flat, square region of a surface of area A ) 4L2 on which have been placed N two-dimensional Gaussian functions of equal width, σ; random (positive and negative) heights, Ci, chosen uniformly in the range -C < Ci < C; and random locations of their centers (xi, yi) chosen uniformly in the ranges -L < xi < L and -L < yi < L. The width is taken to be small compared to the dimensions of the slab, σ , L. The height of the resulting corrugated surface at any position (x, y) is given by -(x - xi)2-(y - yi)2 Ci SN(x, y) ) e 2σ2 σ(2π)1⁄2 i)1 N
∑
(A1)
We define the density of Gaussian functions F ≡ N/4L2 and the density-scaled surface -(x Ci e 1⁄2 σ(2πF) i)1 N
S(x, y) )
∑
xi)2-(y
-
1 2πFσ2
(A2)
2σ2
〈∑
-x i2-y i2 N
N
Cie
2σ2
i)1
∑
-x i2-y i2
Cje
2σ2
j)1
〉 (A3)
Since Ci and Cj are uncorrelated, A3 simplifies to
〈S2 〉 )
1 2πFσ2
〈
N
∑ i)1
(
Ci2 exp -
x i2 + y i2 σ2
)〉
)
S′2x + S′2y ) sec2 R - 1 Define σS2 ) C2/12σ2 and W ) (Sx′/σS)2 - 1)/σS2. The probability distribution of
C2 6
(A4)
using the fact that Ci, xi, and yi are also uncorrelated and can be averaged independently. Invoking the central limit theorem, we see that the height of the surface at each position will be normally distributed with zero mean and variance of C2/6. We denote the slopes of the surface contour, that is, the partial derivatives in the x and y directions, as Sx′ and Sy′, respectively. It is important to note that for sufficiently large N, the gradients are uncorrelated. The mean values of Sx′ and Sy′ are zero. The means of the squares are given by
(A6)
+ (Sy′/σS) ) (sec2 R the sum of squares of two independent standard normal variates is a χ2 distribution with two degrees of freedom. 2
1 W P(W) ) χ22(W) ) e- 2 2 ∂W sin R -(sec2 R-1)⁄2σS2 e P(R) ) P{W(R)} ) ∂R σ 2 cos3 R S
(A7)
Substituting σS2 ) C2/12σ2 and observing that the distribution is uniform in azimuthal angle β for normal incidence, we obtain
Pnormal(R, β) )
yi)2
As shown below, the statistical properties of S(x, y) are independent of N in the limit of large N. We examine the statistical distribution at the position (0, 0) only, as all positions far from the edges are randomly distributed and equivalent. The mean surface height, 〈S(0, 0)〉 ) 〈S〉 ) 0, since the Gaussian heights, Ci, are symmetric about zero. The mean square surface height, 〈S2〉, is
〈S2 〉 ) 〈S2(0, 0) 〉 )
b |-1] ) arccos[(S′2x + S′2y + 1)-1⁄2] R ) arccos[|Q
6σ2 sin R -6σ2(sec2R-1)⁄C2 e πC2 cos3 R
(A8)
The mean value of tan R can be obtained by integrating eq A8:
A ) 〈tan R〉 )
C π σ 24
1⁄2
( )
(A9)
A is the measure of surface roughness employed throughout the paper. Substituting eq A9 into eq A8 gives
Pnormal(R, β) )
1 sin R -π(sec2R-1)⁄4A2 e 4A2 cos3 R
(A10)
Thus, the distribution function is completely characterized by a single parameter, the roughness A. In the case of normal incidence, a group of trajectories incident on an area dx dy projected down onto the global surface plane will be equally likely to hit this area, regardless of the local surface structure. For off-normal incidence, the distribution derived above must be reweighted, as the probability of hitting an area dx dy changes with both the incident angle and the local normal direction. In Figure 1, we consider a line of sight along the x-axis and an incoming trajectory with direction given by ˆ defines the direction of the local unit vector ˆI. The unit vector Q normal at the point of impact. The polar angle of the trajectory, θi, is measured from the global normal, Nˆ. From Figure 1, in 2D, the reweighting factor to account for a tilted local normal or off-normal incidence is
dx cos(θi - R) ) dx cos R When generalized to 3D, eq A11 becomes
(A11)
Accommodation of Gases at Rough Surfaces
ˆ Iˆ · Q ) sin θi tan R cos β + cos θi ˆ · Nˆ Q
J. Phys. Chem. C, Vol. 113, No. 6, 2009 2367
(A12)
where β is the azimuthal angle of the local normal, and the azimuthal angle of the incident trajectory has been defined as φi ≡ 0. Equation A12 can take negative values, which correspond to areas of the surface that slope downhill at an angle greater than the incident angle, so the particle does not strike the surface at this point. Thus, we restrict the weighting factor to positive values, giving the final result for the (un-normalized) probability distribution:
PA(θi ;R, β) )
1 sin R -π(sec2R-1)⁄4A2 e · 4A2 cos3 R Max[{sin θi tan R cos β + cos θi}, 0]
(A13)
This expression is not rigorously correct for two reasons. First, those trajectories that do not strike the surface because the weighting factor of eq A12 is negative will strike the surface at a nearby position. Equation A13 assumes this nearby point will be uncorrelated with the initial impact point, which is not necessarily true. Second, the distribution function of eq A13 assumes that, when shadowing occurs (i.e., when large bumps screen the immediate area beyond from impact by a grazing angle trajectory), the area screened exhibits a statistical distribution of local normals. To test the importance of these simplifying assumptions, we have computed PA(θi; R, β) by numerical simulations of an exhaustive number of trajectories. Even for highly corrugated surfaces and grazing incidence angles, eq A13 reproduces the numerical results to remarkable accuracy. Appendix B: Molecular Dynamics Simulations The molecular dynamics simulations were performed by standard procedures.22 The (111) face of the platinum crystal was represented by a four-layer slab of 256 atoms, with 64 platinum atoms in each layer. The atoms in the top three layers move according to classical mechanics, whereas those in the bottom layer were held fixed in their (0 K) equilibrium positions (Pt-Pt distance ) 2.7746 Å). Periodic boundary conditions were imposed in the x and y directions (z is the surface normal direction.) The platinum atoms were assumed to interact via nearest-neighbor harmonic forces, with a force constant of 262 kJ/mol-Å2. This is sufficient for roughly approximating the average phonon frequency of bulk platinum, but it does not accurately represent the detailed phonon spectrum; neither does it necessarily reproduce the surface atom motions quantitatively. This deficiency is not important for the current purpose of comparison between two model calculations. The Xe-Pt interaction potential was taken to be a pairwise additive sum of Lennard-Jones potentials, with parameters σ ) 5.5 Å and ε ) 0.843 kJ/mol. A cutoff distance of 9.893 Å was employed.
In addition to the Lennard-Jones pairwise terms, an xyindependent “background” potential was included. The background potential was chosen to be of the form
Vbackground ) -A exp(-Bz2)
(B1) -2
where A ) 26.385 kJ/mol and B ) 0.01813 Å . The background potential introduces the well-known “smoothing” of the corrugation at metal surfaces due to the delocalized conduction electrons.23 The combination of Lennard-Jones and background potentials reproduces the measured 26.7 kJ/mol isosteric heat of adsorption of xenon on Pt(111) at zero coverage. Trajectories were integrated with a time step of 1.0 fs. A total of 2000 trajectories were run for each comparison in Figure 2. Trajectories that remained trapped on the surface for longer than 20 ps were assumed to be completely accommodated at the surface temperature. Acknowledgment. This research was supported by the National Science Foundation MRSEC Center for Interface Structures and Phenomena, DMR-0520495; National Science Foundation Grant CHE-0314208; and the Air Force Office of Scientific Research MURI Center for Materials Chemistry in the Space Environment, AFOSR F49620-01-1-0335. References and Notes (1) Goodman, F. O.; Wachman, H. Y. Dynamics of Gas-Surface Scattering; Academic: New York, 1976. (2) Logan, R. M.; Keck, J. C.; Stickney, R. E. In Rarefied Gas Dynamics; Brundin C. L., Ed.; Academic: New York, 1967, p 49. (3) Logan, R. M.; Stickney, R. E. J. Chem. Phys. 1966, 44, 195. (4) Grimmelmann, E. K.; Tully, J. C.; Cardillo, M. J. J. Chem. Phys. 1980, 72, 1039. (5) Tully, J. C. Ann. ReVs. Phys. Chem. 2000, 51, 153. (6) Head-Gordon, M.; Tully, J. C. J. Chem. Phys. 1995, 103, 10137. (7) Tully, J. C. J. Chem. Phys. 1990, 92, 680. (8) Goodman, F. O. J. Phys. Chem. Solids 1965, 26, 85. (9) Trilling, L.; Hurkmans, A. Surf. Sci. 1976, 59, 361. (10) Barker, J. A.; Auerbach, D. J. Chem. Phys. Lett. 1979, 67, 393. (11) Steinbruchel, Ch. Chem. Phys. Lett. 1980, 76, 58. (12) Steinbruchel, Ch. Surf. Sci. 1982, 115, 247. (13) Sitz, G. O.; Kummel, A. C.; Zare, R. N.; Tully, J. C. J. Chem. Phys. 1988, 89, 2572. (14) Berenbak, B. Z.; Riedmuller, B.; Papageorgopoulos, D. C.; Stolte, S.; Kleyn, A. W. Phys. Chem. Chem. Phys. 2002, 4, 68. (15) Carlsson, A. F.; Madix, R. J. Surf. Sci. 2000, 458, 91; 62, 470. (16) Kindt, J. T.; Tully, J. C. Surf. Sci. 2001, 470, 149. (17) Yan, T. Y.; Hase, W. L.; Tully, J. C. J. Chem. Phys. 2004, 120, 1031. (18) Yan, T. Y.; Hase, W. L. J. Phys. Chem. B 2002, 106, 8029. (19) Cook, S. R.; Hoffbauer, M. A.; Clark, D. D.; Cross, J. B. Phys. ReV. E 1998, 58, 492. (20) Cook, S. R.; Hoffbauer, M. A. Phys. ReV. E 1997, 55, R3828. (21) Nathanson, G. M. Annu. ReV. Phys. Chem. 2004, 55, 231. (22) Muhlhausen, C. W.; Williams, L. H.; Tully, J. C. J. Chem. Phys. 1985, 83, 2594. (23) Zaremba, E.; Kohn, W. Phys. ReV. B 1977, 15, 1769.
JP8077634
