Langmuir 1992,8, 2790-2803
2798
Melting Transition of Ethylene on Graphite Ailan Cheng' and Michael L. Klein Department of Chemistry and Laboratory for Research on the Structure of Matter, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 Received March 25,1992. In Final Form: September 11,1992 The melting behavior of submonolayerpatches of physisorbed ethylene molecules on the graphite b a d plane has been investigatedby molecular dynami~calculations.The simulationresultareveal the significant role played by the substrate in determiningthe behaviorof the incommensurateoverlayer. The competition between the ethylene-ethylene intermolecular interactions and the two-dimensional solid-substrate interactionscauses the solid overlayer to be farther away from the surface than in the case of a liquid film. This competition may be related to the observed continuous melting transition. In the transition region, the adsorbate forms solid patches connected by liquidlike domain walls. The coherence (correlation) length of the ordered domains and the structural relaxation time are estimated.
Introduction
1 / I -oln
coexistence
In spite of intensive study, the melting transition in two dimensions is still a controversial problem. Numerous experiments have been carried out to investigate this phenomenon. Many systems,l such as Kr, Nz, and CH4, exhibit fmt-order meltingtransitions in the submonolayer region; however, under similar conditions Ar and CZH4 melt continuously. Two-dimensionalmelting has attracted many theoretical studies. The KTNHY theory, proposed by Kosterlitz and Thouless,2 later modified by Nelson and Halperin? and independently by Young? predicted that melting is a two-step process involving the unbinding of two types of topological defects: the dislocation and disclination. Thus, the two-dimensional solid is predicted to melt first to an anisotropic liquid, and then to an isotropic liquid. Neutron diffraction6t6has been used to compare the melting behavior of CH4 and CzK on graphite. Molecular dynamics simulations have been formed to investigate the melting behavior of methane on graphite.' Here, we concentrate on the melting transition of ethylene physisorbed on graphite. In the submonolayer region ethylene forms an orientationally ordered solid phase (OLD)below 40 K (see the phase diagram in Figure 1). The unit cell, which is rectangular, contains two molecules with the C = C bond of the molecule parallel to the surface. As temperature increases an orientational disordered phase (DLD) is observed, in which the C = C bond begins to rotate about an axis normal to the surface and the molecules form a triangular lattice. This triangular phase melts continuously at 68 K.&" It is an intriguing question as to whether or not the melting transition of CZH4 on graphite is of the KTNHY type. The latest low-energyelectron diffraction (1) Strandbug, K. J. Reo. Mod. Phys. 1988,60, 161 and references therein. (2)Kosterlitz, J. M.; Thoula, D. J. In Progrees in Low Temperature Physics; Brewer, D. F.,Ed.; North-Holland: Amsterdam, 1973;p 373. (3)Nelaon,D. R.; Halperin,B.I. Phys. Reu.E 1979,19,2457.Halperin, B. I.; Nelson, D. R. Phye. Reo. Lett. 1978,41,121. (4)Young, A. P. Phye. Reu. E 1979,19,1856. (5)h e m , J. 2.; Paeeell,L.; Heidemann, A. D.; Richter, D.; Wicknted, J. P. Phys. Reu. Lett. 1988,61,432. (6)Satija, S. K.; Paasell, L.; Eckert, J.; Elleneon, W.; Patterson, H. Phys. Reo. Lett. 1983,51,411. (7)Kim, H.-Y.; Steele, W. A. Phys. Reo. E 1992,45,6226. (8) Zhang, Q.M.; Feng,Y. P.; Kim, H. K.; Chan, M. H. W. Phys. Reo. Lett. 1986,57,1456. (9) Kim, H.K.; Zhang, Q. M.; Chan,M. H. W. Phys. Reo. Lett. 1986, 56, 1579. (10)Larew, J. 2.; Paasell, L.; Ravel, B. Can. J . Chem. 1988,66,633. (11) Larese, J. 2.;Rollefson, R. J. Surf. Sci. 1983,127,L172.
0743-7463192/2408-2798$03.oOlO
h
d
1.0
v
-
e,
2 09 L g 08
01-0*
e,
OCD
[
OH) *DLD
/
u
07
06
20
30
40
50
60
70
80
'T(K)
Figure 1. Phase diagram of ethylene on graphitetaken from ref 12. Symbolshave the following meanings: LD, low density;HD, high density; 0, ordered; D, disordered; C, commensurate; I, incommensurate; V, vapor; and L, liquid.
(LEED) experiment12indicated that there exists a commensurate orientationally disordered phase (CDLD) between the triangular and fluid phases. Computer simulation techniques have previously been used to study ethylene physisorbed on graphite. "6 and Klein13 used a constant-pressure two-dimensional version of the Parrinello-Rahman algorithm to inveetigate both the low- and high-density phases of C2)4 on a flat surface for a range of temperatures. They predicted the existence of the OLD phase, later observed in the neutron scattering.1° Subsequent work" included the effect of substrate corrugation. These simulationresultd4seemed to confirm the melting behavior observed in experiment, but the calculated melting transition temperature was much lower than the experimental one. In an effort to better understand the continuousmelting behavior and to obtain a realistic transition temperature, we have carried out new molecular dynamics simulations. We employ slightly different potential models. In particular, a more realistic anisotropic potential's is used for the adsorbatesubstrate interaction. The latter r e f l a the fact that carbon atoms in the graphite are more polarizable in the basal plane than in the direction perpendicular to it. Effectively, this contribution to the surface corrugation almost doubles the value of the ueual (12)Eden, V. L.;Fain,S. C., Jr. Phye. Reu. E lSS1,43,10697. Eden, V. L.; Fain,S. C., Jr. Preprint, 1991. (13)No&, S.;Klein, M. L. Phys. Reo. Lett. 1984,63,818. (14)Moller, M. A.;Klein, M. L. Chem. Phys. 1989,129,235. (16)Vidali, G.;Cole, M. W. Phy8. Reo. E 1984,29,6738.
Q 1992 American Chemical Society
Langmuir, Vol. 8, No.11,1992 2799
Melting Transition of Ethylene on Graphite Table I. Atom-Atom Lennard-Jones Potential Parameters for the Adsorbate-Substrate Interactions Udr) = 4e((u/r)'f - (u/r)9 i-j
B
C-C (substrate)
H-C (substrate)
(K)
47.72 17.00
(A)
C-C
3.30 2.98
C-H
a
Table 11. Atom-Atom (exp-6) Intermolecular Potential Parameters (Model I) Odr) = A exp(-ur) - B/Z i-j A (kJ/mol) B (kJ/mol) Q (A-9 ~~
C-C C-H
H-H
311 540.6 39 375.6 16 736.0
-2238.4 -581.6 -150.6
Table 111. Atom-Atom (exp-6) Intermolecular Potential Parameters (Model 11). i-j A (kJ/mol) B (kJ/mol) u (A-1)
3.60 3.67 3.74
isotropic potential,16J7used in the previous study.14 In the simulation of N2 on graphite,l8themelting temperature was raised by 10 K when the anisotropic substrateadsorbate interaction was employed. We also explicitly include the intermolecular quadrupole-quadrupole interaction in a second set of simulations. Details of the Simulations In order to avoid undesirable constraints on the system due to the periodic boundary condition applied in the simulation and to create an adsorbed system close to the experimental conditions, a patch with a free boundary has been studied here. A large rectangular simulation box, 48a X 26(3)'l2a, where a = 2.46 A is the lattice parameter of the substrate, has been used to reduce edge effects. The patch contains 480 molecules and is about 100 A X 100 A. This is comparable with the reported experimental domain size of about 200 A5vg Each ethylene molecule is taken as rigid, with a C=C bond length of 1.337 A, C-H bond length of 1.086 A, and H-C-H angle of 117.4'. The six atoms in eachmolecule interact with those in other molecules and the substrate via pairwise additive potentials. The interaction of the C2H4 molecules with the substrate is considered to be of the Lennard-Jones (12-6)type. The substrate interaction also included the anisotropiceffect mentioned above, with parameters YA = 0.4 and YR = -0.9.15J8 The adatomsubstrate potential parameters were taken from the previous work14 and are listed in Table I. With this potential,the amplitudesof the corrugationfor the C atomsubstrate and H atom-substrate potentials are 56 and 6 K,respectively; these values are almost twice as large as those for the isotropic potential. Here, the corrugation is defined as the differencebetween the potential maximum and minimum across the surface at the minimum of the holding potential. The corrugationexperienced by a C2H4 molecule varies with the orientation of the admolecule and ranges from 10 to 130 K. The intermolecular potential is built from interactions of pairs of atoms of the exp-6 form. The first set of parameters used (model I) was taken from a previous study;14they are listed in Table 11. The second, model 11, also included electrostatic terms, and was taken from the work of william^.'^ The nonbonded van der Waals interaction (listed in Table 111)is different from that of model I. The molecular charge distribution is modeled by charges located on each of the six atoms (QH = 0.153 e, qc = -0.306 e). This distribution gives Qxx = -2.84 B, Q,, = 0.95 B, and Qzz = 1.89 B, when the molecular fixed principal axes are chosen the following way: z axis is (16)Stele, W.A. Surf. Sci. 1973,36,317. (17)Carlos, W.E.;Cole, M. W. Surf. Sci. 1980,91, 339. (18)Joshi, Y. P.;Tildesley, D. J. Mol. Phys. 1985,55, 999. (19)Williams, D.E.,Starr, T. L. Comput. Chem. 1977,1,173.
H-H
367 250.6 65 485.0 11 677.0
-2414.0 -573.0 -136.0
3.60 3.67 3.74
pc = -0.306 e and p~ = 0.153 e were used as atomic charges.
parallel to the C = C bond, x axis is perpendicular to the molecular plane, and they axis is parallel to the molecular vary somewhat with plane. Experimental values for QL2@ the method of measurement. However, collision-induced absorption data suggest Qxz = -3.25 B, Qyy = 1.62 B, and in fair accord with our model. Q Z E = 1.63 The molecular dynamics calculations were carried out using the conventional microcanonical method. The equations of motion for the translational degrees of freedom were solved by a fifth-order Gear predictorcorrector algorithm.21 A fourth-order algorithm was used for the rotational motion. The coordinates were updated every 2 fs. The simulationstarted from a triangular lattice with the C = C parallel to the surface and herringbone orientations. Typically,the systems were first allowed to relax for 20 ps or longer. Then, configurations were collected over longer runs for subsequent analysis. In the region of the melting point, run lengths greater than 200 pa were required to ensure that the system was equilibrated. Results of Calculations A. Model I. Simulations were carried out for several temperatures. Thermodynamic data and information on the orientationalorder were obtained from the trajectories. Analysis of configurations for different temperatures indicates an increase of translational disorder as the also temperature rises. The total potential energy, Umd, increases with temperature. The change in slope between 40 and 50 K, seen in Figure 2, is likely the melting transition. Diffusion coefficents for translational motion of the center of mass and rotational motion about the three principal axes have been calculated from the velocity and angular velocity correlation functions.22The orientation of each molecule relative to the substrate is described by the Euler angles: 8,the polar angle,4, the azimuthal angle, and $, the angle of rotation about the C = C axis of the molecule. Here, the vector normal to the substrate is chosen to be the L axis of the space-fixed frame. Translational motion is separated into motion parallel to and perpendicular to the surface. Only the diffusion coefficient for motion parallel to the surface is shown in Figure 2. T h e motion perpendicular to the surface (z direction) is almost harmonic-oscillator-like;there is no layer promotion, and therefore there is no diffusion along the z direction. It is not our purpose here to stress, and it is impossible for us to identify, the order of the melting transition since the simulations are carried out on such a small system (N= 480) and span a rather short (less than nanosecond) time scale. Instead, we concentrate on grow changes of the order as the temperature rises. The logarithm of the translational diffusion coefficients can be fitted by two (20)Dagg,I. R.;Read, L. A. A.; Smith, W. Can.J. Phys. 1982,60,1431. (21)Gear, C. W.T h e Numerical Integration of Ordinary Dfierential Equations of Various Orders. Report ANL 7126; Argonne National Laboratory: Argonne, IL, 1968, Numerical Initial Value Problem in Ordinary Differential Equations; Prentice-HalL Eaglewood Cliffa, NJ, 1971. (22)Allen, M.P.;Tildealey, D. J. Computer Simulation of L i q u i l , Clarendon: Oxford, 1987.
Chew and Klein
2800 Langmuir, Vol. 8, No. 11, 1992
0
-
h
-24
: X v
1-25
0
0 0
-
3
-26
h
I
0 0
2
YI
"E
0
P
21 -
v
h
n
0
0
I
I
.
0 0
-
5P 1
0
-
0 0 0 A
n
2
A A A
=,
, I
0
I
60
45
75
T(K)
Figure 2. Simulation results based on the intermolecular potential model I without electrostaticinteractions. The upper panel contains the temperature dependenceof the configurational energy. A change of slope occurs around 47 K. Diffusion coefficients are shown for translation (centerpanel) and rotation (lowerpanel). Triangles,squares, and dotsarefor rotations about the three principal molecular axes x (normal to the molecular plane), y, and z (along the C=C bond), respectively. In this case, the x axes of most molecules are parallel to the surface normal. straight lines, corresponding to the low- and high-temperature regions, respectively. The activation energies (Ec)obtained for the two regions differ significantly: E, = 240 K for the low-temperature region and E, = 65 K for the high-temperatureregion. We identifythe temperature at which there is a change of the adsorbate mobility, 47 f 5 K, as the melting point. As in previous work, this temperature is much lower than the experimental value (68 K). Although one cannot really compare directly with the previoussimulation,14which was carried out at a higher coverage than the current system, it seems there is no significant improvement over the early work even when an anisotropic surface potential is included. In the case of monolayer N2 on graphite which forms a d 3 X d3 commensurate phase at low temperature,18the anisotropic potential (effectivelya higher substrate barrier) naturally stabilizesthe solid phase and therefore raises the transition temperature. However, in the present case, CzHd forms an incommensurate phase and the effect of the substrate corrugation is averaged out. Evidently, the higher substrate barrier used here does not stabilize the solid phase as in the case of the commensurate Ndgraphite system. Below the melting temperature, the C2H.4molecules form a triangular lattice as indicated by the pair correlation function for the centers of mass. The intermolecular spacing is 5.0 A. The nearest-neighbor distance is much larger than that obtained from neutron experiments (4.65 A1.596 In the simulation there is no noticeable change of the nearest-neighbor distance with temperature. The
1
30
.=
A: 45
D
60 T (to
Figure 3. Simulation results based on model I1 that include electrostatic interactions. "he upper panel containe the temperature dependence of the configurationalenergy. A change of slope occurs around 57 K. Diffusion coefficients are shown for translation (centerpanel) and rotation (lower panel). Triangles, squares, and dots are for rotations about three principal axm x , y, and z, respectively. distribution of the three Euler anglesrevealsthat the c--C bond always lies parallel to the surface. The angle distribution, p($), indicates that the molecular plane stays parallel to the surface. The azimuthal angle, 4, shows that there are two preferred orientationswhich correspond to the herringbone structure, though the molecules can rotate about an axis normal to the surface as well. The angle between the two sublattices is 60°. B. Model 11. In the studies of N2 on graphite, the quadrupole moment plays an important role. Ethylene has a much larger quadrupole moment than the N2; a natural first step to improve the potential model ia to include the quadrupole-quadrupole interaction. This contribution is included here by fractionalcharges located on each atom as described above. The calculations were repeated with the potential parameters given in Table 111. When the quadrupolequadrupole interaction is included in the simulation, the C2H.4X2H.4pair interaction becomes stronger in comparison with model I. Typically, the electrostatic terms contribute 10% of the CzH&zH4 pair interaction. The translational order increases significantly, compared with the corresponding data without the electrostatic interactions. Diffusion coefficients are shown in Figure 3. However, there is no apparent difference in rotational diffusion coefficients. The simulation results indicate that disorder develops gradually. Changes of potential energy with temperature suggest that melting occurs below 60 K (Figure 3). The transition temperature deduced from fitting diffusion Coefficients to an Arrhenius law is 57 K. The activated
Melting Transition of Ethylene on Graphite
Langmuir, Vol. 8, No. 11, 1992 2001
X
v
-19.4
I'
0
,I'
3
4
z(R)
Figure 4. (Upper panel)Temperature variation of the ethylene substrate interaction energy, Up,based on model 11. The anomalous temperature dependence is due to the fact that the solid overlayer is farther away from the surface than the liquid film. (Lower panel) Center of mass density distribution along the surface normal (z direction).
energies obtained are E, = 400 K and E, = 200 K for lowand high-temperature regions, respectively. The stronger intermolecular interaction has raised the melting point. The low-temperature structure is a triangular lattice with a smaller intermolecular separation (4.7 A) than that obtained in model I, due to the extra attractive electrostatic interaction. This new value is in much better agreement with experiment (4.65 The distribution of the polar angle and azimuthal angle did not change much. In the low-temperature phase, molecules tilt 40° away from the surface to minimize the adsorbate-adsorbate interaction, especially the quadrupole-quadrupole interaction. Note the quadrupole-quadrupole interactions between linear NZmolecules prefer herringbone orientations. Energetically, however, the substrate favors the configurationsin which molecules are parallel to the surface. As the temperature increases, lattice expansion weakens the quadrupole-quadrupole interaction and the orientation of the molecular plane changes from a tilted position to one parallel with respect to the surface. The diffusion coefficient for rotation about the C==Caxis is greater than the other components; however, the azimuthal angle distribution shows that molecules have no preferred orientations a t 60 K, while there is clear preference for the molecular planes to be parallel to the surface, even at 70
K. The changes of average molecular orientation result in an interesting temperature dependence of the C&substrate interaction energy shown in Figure 4. The intermolecular interactions become weaker as the temperature rises, but the CzH4-graphit.einteraction becomes stronger. The turning point coincides with the melting point. A detailed analysisof the density distribution along the surfacenormal, p ( z ) , confirma that the solid layer stays farther from the surface than the liquid layer. The competition between the CzH4-C& interactionand solidsubstrate versus liquid-substrate interactions influences the melting transition. In neutron scattering, a solidlike contribution was observed above the melting transition temperature.= There is substantial translational diffusion at temperatures (23)Larese, J. 2.; Passell, L. Privata communication, 1991.
Figure 5. Trajectory plots for the molecular center of mam covering 16 ps based on model 11at T = 40,50, and 55 K for the top to bottom panels, respectively. In the latter, the system
consists of several ordered domains.
well below the melting point in the simulation results. Close examination of configurations (Figure 5 ) near the melting region suggests that the adsorbate forms small solid patches, with liquidlike domain walls. The timedependent configurationsdisplayedon a graphics terminal show that these domains move relative to each other and some molecules stay close to their initial position for a certain period of time, and then move away. Therefore, we monitored the number of molecules that stay close to their lattice sites as a function of time, and this is shown in Figure 6. The number of molecules which have a rootmean-square displacement smaller than 15% of the intermolecular distance decays exponentially with time. The short-lived solid component is somewhat similar to what was observed in neutron scattering. As temperature increases, the time-dependent trajectories show that the solid domains shrink and liquidlike domain walls grow; eventually the system melts. The average size of the patches, or correlation length, can be
C h e w and Klein
2802 Langmuir, Vol. 8, No. 11, 1992 500
I
I
300 h
Y 4
i
200
100
0
10
0
30
20
40
NPS)
Figure 6. Number of molecules which have root-mean-square displacements less than 15% of the nearest-neighbor distance for simulations based on model 11. The curves from top to bottom are for T = 40,51, 53, and 60 K,respectively. The inset shows the relaxation times obtained from the N&) plots.
L 0
5
1
I
I
10
15
20
25
r(4 Figure 7. Bond correlation function g&) based on model 11, defined in the text, for T = 51,53, and 60 K for the top to bottom curves, respectively. The inset shows the correlation lengths for the correspondingtemperatures. The lowest temperaturesample is fully ordered.
obtained from g&), the bond correlation function,defined -24
,SiCk+'jd
gs(r) =
)
where m,k and j, 1 denote pairs of atoms separated by a distance r and t9,k is the angle made by the vector from atom m to atom k, relative to some arbitrary coordinate axis. If the lattice has perfect hexagonal symmetry, g d r ) will be unity at each r correspondingto a peak in g(r).We have evaluated g d r ) at several temperatures. The longrange decay behavior, fitted to an exponential function, exp(-rlU, yields the correlation length [. As shown in Figure 7, the correlation length in the low-temperature (40 K)solid phase is comparable with the system size (150 A in diameter), and agrees well with experiments. The correlationlength becomes smaller as the lattice melts. At 60 K, only the nearest-neighbor molecules are correlated. ~~
~
~~
The correlation length obtained in the neutron experiment changed from 200 to 30 A between 50 and 70 K.6J3 In the previous molecular dynamics (MD) study of higher density systems (n = 0.8 and 0.91, the melting transition seemed to be induced by the moleculesstanding up.14 The current simulation results did not shown substantial numbers of erect molecules. But their number increases when the temperature reaches 50 K. The previous simulation14was carried out at a much higher coverage and a full monolayer. In that caee, thermal expansion can only be realized when molecules stand up to generatethe space or by layer promotion. In the current simulation, only 70% of the substrate is covered. Even at 120 K,70% of molecules stay parallel to the surface. Standing up is not energetically favorable, due to the substrate attraction. The single molecule-substrate interaction varies significantly with the orientation of the admolecule. The lowest energies are -2600, -2000, and -1710 K, respectively, for the three orientations lying parallel to the surface, standing on a side (edge), and standing on end. In the recent study of rod-shaped molecule~,2~ it is found that tilting is a more energetically favorable mechanism for vacancy formation in the butane monolayer. From the diffusion coefficient data, it seems that translational disordering is associated with change in orientations of the molecules. The orientational distributions of the Euler angles are calculated for the solidlike and liquidlike components separately at 50 K (Figure 8). The distribution of the /3 angle is almost the same, indicatingthat the erect molecules are not the major cause of the disorder. Projections of the molecular Cd: bond on the x-y plane show some preferred angles when molecules are in the solidlike region; however, molecules in the liquid region can rotate much more freely about the axis normal to the substrate. Moleculesin the solid region tilt 40° awayfrom the surface;the width of the distribution is about 20°. The angle distribution p(#) for the liquid component has wider peaks centered at fWO, which correspond to the molecular plane being parallel to the substrate. The liquid part exhibits more orientational disorder than the solid part. As temperature increases, the in-plane rotation and the libration of molecules about their C==C axis become much easier. Some molecules can stand up and hence, temporarily, create more room for the nearby molecules. Therefore, there is higher local mobility when moleculesstand up. Effectively,these local defects may induce melting. For the purpose of analysis we again separate molecules into two groups, those in the solidlike region and those in the liquidlike region. The average peak position in the density distribution p(z) for the solid part is 0.1 A higher than that of the liquid component at 50 K. The average solid-substrate interaction is 57 K weaker than that between the liquid and substrate.
Conclusion New molecular dynamicssimulationshave been carried out for submonolayer films of ethylene on graphite. The results indicate that the melting transition occurs gradually, which is in good accord with the experimental fih1dings.~*~~~9 In the transition region, the adsorbate forms solid domains surrounded by liquidlike walls. Molecules in the solidlike regionsstay closeto their originalpositions for a certain time, and then move away, to joint new patches. This kind of transient solid domain may explain
~~~~
(24) McTague, J. P.; Frenkel, D.; Allen, M. P. In Ordering in Two Dimensions; Sinha, Ed.; North-Holland Amsterdam, 1980; p 147.
(26) Hausen, F. L., Taub, H. Phys. Rev. Lett. 1992,69, 662.
Melting Transition of Ethylene on Graphite .2
I
I
Langmuir, Vol. 8, No. 11, 1992 2003
I
I
I
10
.03
h
Q
.02
.01
0
I
I
-100
I
0 #(del31
I
I
100
.03
*
h
.04 p(molecule/A')
.02
.06
Figure 9. Isosteric heat of adsorption from the experimental study (dots)and current simulation based on model I1 (quares) at 183 K. The C Z H ~ - C Zinteraction H~ in the simulation seems to be tooweak and the ethylene-substrateinteractiontoo strong, compared with experiment. Table IV. Molecular Dynamics Results for the DLD Phaw, Baaed on the Two Intermolecular Potentials lattice constant (A) T, (K) tilt anale0 (den) model I model I1 experiment
5.0 47 4.7 51 4.65 68 a Tilt of molecular plane away from the surface.
.02
Y
P
-
0
8 v
I
I
.01
0 -100
0
100
l#(deg)
Figure 8. Probability distributions p(/3), p(#), and p(+) based on model II at T = 60 K. The curves are for solid (solid lies) and liquid (dotted lines)components averaged over 10 ps. The distribution p(B) indicatesthat the molecular C = C axis in both regions lies parallel to the surface. Molecules in the solidlike ) and regions are more ordered in the x-y plane ( ~ ( 4distribution) tilt 40" away from the surface (p(+) distribution). In the liquid regions, molecules rotate almost freely about an axis normal to the surface and molecular planes librate about a plane parallel to the surface. the elastic component seen in a recent neutron scattering study.29 The estimated spatial correlation length of the domains agrees well with the experimental results. The analysis of the orientational motions indicated that most molecules have their C==C axes parallel to the surface in the submonolayerregion. There is no substantial number of molecules standing up. Also, the phenomenon of layer promotion is not observed due to the strong substrate attraction. The rotations about an axis normal to the surface and about the C - C bond increase rapidly in the transition region. The tilting of molecules in the solid lattice c a w the layer to stay farther away from the surface than molecules in the liquid layers which lie parallel to the surface. The substrate interaction therefore favors the liquid phase, energetically. This effect causes an anomaloustemperature dependenceof the Cz&-substrate interaction, which could perhaps be related to the continuous melting behavior in adsorbed ethylene. When electrostaticquadrupole-quadrupoleinteractions are included in the simulation, we obtain structure and orientational order in better agreement with the experimental The melting transition temperature is raised by 10 K, but it is still lower than the experimental value. The lower transition temperature is, in part, due to the effect of the small sample size used in this study.
0 40 35
The comparison between results of two potentials (i.e., without and with electrostatic interaction) are listed in Table IV. To further test the potential models, we have estimated the isostericheat of adsorptionz6at 183 K to compare with the experimental data at the same temperat~re.~'Both the experimental and simulation values are plotted in Figure 9. A natural conclusion is that the C2H4-Czfi interaction is still not strong enough. One can imagine a further increase in either the quadrupole moment or the van der Waals interaction and at the same time modification of the repulsive core to produce the correct lattice constants. In fact, a radically different intermolecular potential model, derived from the ab initio calculations, was used in the study of dynamical properties of solid ethylene.z8 For this model, the point charges are shifted away from atomic sites to obtain a better representation of the calculated electrostatic interactions between molecules. In addition the C-C interaction well depth is 20 K deeper than the value used in the current study. Perhaps, by using this potential, one can raise the melting temperature, but the potential effectively has 12 interaction sites which will increase the computation time dramatically. Clearly, much still remains to be done to obtain a quantitative modeling of the ethylene/graphite system. Acknowledgment. We thank Sam Fain, John Z. Larese, and Larry Passell for stimulating us to initiate this work, for sharing their unpublished results, and for valuable discussions. This research was supported in part by the NSF and the donors of the Petroleum Research Fund, administered by the American Chemical Society. (26) Steele, W. A. The Interaction of C u e s with Solid Surfaces; Pergamon: New York, 1974. (27) Avgul, N. N.; Kiselev, A. V. Chem. Phys. Carbon 1970,6,1. (28) Luty, T.; van der Avoird, A.; Berns, R. M.;Waeiutynski, T. J. Chem. PhYS. 1981, 75,1451.
