- Email: [email protected]
The continuous melting transition of ethylene on graphite: A molecular dynamics study
Chemical Physics 129 (1989) North-Holland, Amsterdam
235-239
THE CONTINUOUS MELTING TRANSITION A MOLECULAR DYNAMICS STUDY
OF ETHYLENE
ON GRAPHITE:
Michael A. MOLLER Department of Chemistry, McMaster University, Hamilton, ‘Ontario. Canada LBS 4Ml
and Michael L. KLEIN of Chemistry,University ofPennsylvania, Philadelphia, PA 19104-6323, USA
Department Received
1 May 1988
The melting of monolayer films of ethylene on graphite has been investigated by molecular dynamics calculations. Our results suggest that melting is influenced by the presence of individual molecules standing erect, with their C=C axes perpendicular to the graphite surface. These molecules lower the effective density in their vicinity and because their number increases with temperature, cause the melting to appear continuous.
1. Introduction
graphite identified translation-rotation coupling as playing an important role in the dynamics of ethylene overlayers [ 41. However, the process in which molecules become erect does not appear to be related to such cooperative (collective) phenomena or a phase separation. The fluctuations seem to occur randomly in the sample, and persist for times of order lo-20 ps. We suggest that the observed anomalous melting transition for ehtylene on graphite is intimately associated with this novel process. We now describe the MD calculations which lead us to this conclusion.
The diversity of the melting characteristics exhibited by quasi two-dimensional molecular overlayers, physisorbed on graphite, is quite surprising [ l-3 1. For example, submonolayer films of the simplest hydrocarbon, methane (CH,), show a sharp first-order transition; the heat capacity peak associated with melting is narrow and exceeds 100 k, in height [ 11. The behavior of the aspherical molecule ehtylene (C,H,) is quite different. Here, melting extends over several degrees and the heat capacity peak has a maximum height of only about 2 k, [ 21. The latter observation has been interpreted as evidence for a continuous melting transition [ 2 1. In order to seek a possible microscopic explanation of this observation, we therefore embarked on a molecular dynamics (MD) study of the ethylene-graphite system. An analysis of our MD calculations leads us to conclude that the melting of ethylene films physisorbed on graphite is strongly influenced by molecular reorientation. The local density fluctuations required for melting are enhanced when some of the molecules stand up, with their C=C axes perpendicular to the graphite basal plane. A previous MD study of C,H,/ 0301-0104/89/$ ( North-Holland
03.50 0 Elsevier Science Publishers Physics Publishing Division )
2. Calculations We have employed the same model as used previously to discuss the structure and dynamics of the solid phases formed by ethylene on graphite at low temperatures [ 4,5 1. The intermolecular potential for the ethylene molecules was the sum of atom-atom interactions of the exp-6 type (see ref. [ 61, “the potential set IV”). The molecule-surface potential was derived from atom-atom interactions between the atoms of the ethylene molecules and the atoms comB.V.
236
M.A. Moller. M.L. Klein /The melting ofethylene on graphite
prising the graphite basal plane. In order to include the corrugation, we have followed the method of Steele [ 7 1. Here, .the atom-surface potentials are developed in a Fourier series expansion. The first term, V,,(z), yields the flat surface approximation used in our earlier MD studies of this system [4]. The second term in the Fourier series, V, (x,y,z), contains the effect of the surface corrugation. The variation of V,, as an ethylene molecule moves across the graphite surface is about 50 K for the potentials used here. The optimal configurational energy for a single molecule interacting with the graphite basal plane agrees to within a few percent of the value derived from the measured isosteric heat of adsorption extrapolated to zero coverage, namely - 19.9 kJ/mol [ 51. Inclusion of the corrugation term V, imposes an important constraint on the system to be studied. Since the simulation box had periodic boundaries parallel to the (x,y) plane, the potential had to be continuous at these boundaries. To ensure that this was so, the boundaries were kept fixed, with an integral number of graphite hexagons inside the MD box. Simulations were performed in the (N,E, V) ensemble at two different coverages. We define a coverage of y1= 1 to be a commensurate &X v1? R30” monolayer. One set of calculations employed a box with sides 13$a x 23a, where a= 2.4 A is the lattice constant for the graphite basal plane hexagons. There were 598 hexagons within this box. We used N= 168 ethylene molecules, so the coverage in this case was n = 0.84. This value corresponds to the coverage of a complete monolayer of the so-called low-density phase at low temperature. The remaining simulations were carried out with N= 132 molecules in a box with edges 11 ,/?a x 20a. Since this box contained 440 hexagons, the coverage was n = 0.90. This was the lowtemperature equilibrium coverage for our potential model as determined in previous calculations [ 41. By running simulations at this coverage, we eliminate the possibility that the low melting temperature at n=0.84 and the resulting enhanced motion of the molecules is due to a coverage less than the equilibrium value for our model. Experimentally, at low temperature and at coverages less than 0.84, solid ethylene overlayers adopt an incommensurate structure, expanded with respect to the registered phase, in which the C=C molecular axes are predominantly parallel to the graphite basal plane
[4,8]. Moreover, before it melts, this so-called lowdensity solid transforms (at Tz 35 K) to a rotator phase [ 2,4,8,9]. In view of these facts, our MD calculations were initiated from ordered triangular herringbone structures with all C=C axes parallel to the graphite surface. Since we are interested in melting, systems were prepared at temperatures ranging from 40 to 100 K [ 8 1. In all cases, the ethylene molecules became orientationally disordered in the plane of the surface. A typical simulation, run for about 200 ps using a time step of 2.5 fs, took about 20 h of CPU time on an IBM 3090; it was therefore impractical to investigate a larger system size. To locate the melting temperature, the structure amplitude
S(k)+
cos(k.r))
was evaluated as a running average throughout each simulation. Here, k is one of the three smallest reciprocal lattice vectors of the triangular lattice appropriate to the surface coverage of interest, Yis the distance between the centers of mass of two molecules and the summation is over all such pairs of molecules. The structure amplitude was calculated for each molecule in the system as the origin, and then averaged. The average was also taken over all pairs of molecules and time steps. As defined, S(k) = 1 for a perfectly static triangular lattice and falls to near zero for a translationally disordered liquid. In a solid at temperatures above zero, S(k) is reduced to a value less than unity by the Debye-Waller effect (thermal motions), and static disorder (dislocations). In the following discussion and figures, we have taken the structure amplitude to be the average of all three S(k) values where k=(l,O), (0, l)and(l, 1). We also studied in detail the diffusional motions of the ethylene molecules on the surface. The translational diffusion coefficient, D,, was obtained from the time dependence of the mean square displacements of the molecular centers of mass in the (x, Ji) plane. Rotational diffusion coefficients, D,,, were obtained from the appropriate angular velocity autocorrelation functions. We followed rotation about axes in the molecular plane both parallel and perpendicular to the C=C axis as well as rotation about an axis perpendicular to the molecular plane. Rota-
M.A. Moller, M.L. Klein /The melting of ethylene on graphite
tional motion about C=C is by far the most important. Other quantities monitored included the average tilt of the C=C axis away from the plane of the surface, and the configurational energy, which was separated into a molecule-molecule part E( m-m) and a molecule-surface part E( m-s).
0.6 0.5
ti! 3
0.4
6
0.3
E =
0.2
I I
_I
3. Results
0.0
%
0.6
5
O*
E
0.2 0.0 0.6
2
04 .
5 0.2
0
50
toot5o200250
TM @> Fig. 1. Time evolution of the structure amplitude for ethylene overlayers physisorbed on graphite. The upper curve (a) is for T= 53 K and the lower curve (b) is for T= 63 K. Both calculations were carried out at a surface coverage of n=0.90 (see text).
I
l
I
vO.64
. n=O.90
t
,I
*fl t,,
I
’
b
I
, ’
, * I .
- 04 0.15
t
t:
t
is ii?
8
II
t7; 0.x)
&
I
t (4
9 B 3
8
I’t
0.1 _
Fig. 1 shows the time evolution of the structure amplitude S(k) for points on either side of the melting transition in the system with n=0.90. Close to the melting transition, long simulation trajectories (200 ps) are required to establish equilibrium. In fig. 2a we show the temperature dependence of S( k) for both coverages. The value of S(k) was calculated every 25 steps, and these values were averaged to produce the figure. The error bars represent the standard deviation of S(k). For the series with n=0.84 (dots) there is a decrease in S(k) in the region 5260 K whereas the drop occurs between 57 and 63 K for the series with n = 0.90 (squares). This drop-off
231
t
0.05
t tt + +*
tt I
0.00 35
55
75
95
TEw09 Fig. 2. (a) The temperature dependence of the structure amplitude and (b) the fraction of molecules in the system standing up.
in S(k) identified the melting region for a given coverage. The fluctuations increased considerably as the temperature approached the melting point, representing the formation and annealing of dislocations. Unfortunately, the calculated melting regions are somewhat lower than the experimental counterparts which are centered around 70 and 80 K. This disagreement likely signals inadequacies in the relevant interaction potentials we have employed. Fig. 3 shows the temperature variation of the diffusion coefficients D,, and D,,. A line drawn through the values of D*, extrapolates to zero at about 50 K for n~0.84 and at about 55 Kfor n~0.90. These values lie below the respective melting regions and suggest that D,, is continuous across the melting “transition”. Recent quasi-elastic neutron scattering data are in agreement with this result but also suggest that the melting is associated with enhanced molecular rotation about the molecular axes perpendicular to C=C [ lo]. This latter observation seems to be at variance with our MD results (fig. 3b). This discrepancy may be related to the fact that, in the experiment, the overlayer can expand on heating, whereas the calculations are carried out at fixed coverage, which would restrict the motions of the molecules.
238
M.A. Moller, M.L. Klein /The melting qfethplene
c
I
’
r
I
I
on graphite
-6.2
(a) . n=O.84
w,
t
. .
(al
. n=O.90
-6.9
-7.6
-y rn
9.0
Fg
6.0
0
P
t
t
+k
t
tt
t+:
3.0
+ t t I
0.0 35
A
Ati
b
I
55
T
(;
@*)! @*I 4
I
95
Fig. 3. (a) Temperature variation of the translational diffusion coefficient, D,,. (b)Temperature variation of the rotational diffusion coefficients, D,,,, for nE0.84. Symbols indicate rotation about the C-C axis (circles), and about orthogonal axes perpendicular (triangles), and parallel (crosses) to the molecular plane.
Fig. 4 shows the temperature variation of the configurational energies E(m-m) and E(m-s). These quantities exhibit only a small variation across the melting region. Indeed, crude estimates of the derivatives of these curves with respect to temperature suggest a specific heat anomaly of order 2-3 k, in the relevant transition regions. This behavior is entirely consistent with the measured heat capacity data [ 21. Even though the calculated melting temperature is somewhat low, there is apparent success in modeling the thermal anomaly at the melting transition [ 21. Accordingly, we were encouraged to search for a possible microscopic explanation of this weak transition. Of all the quantities we examined, only one showed any dramatic change at melting and that was the number of molecules standing up in the system. The orientational distribution function for C=C with respect to the surface plane was bimodal, with peaks at 0” and 80” and a well-defined minimum extending from about 20” to 65”. Based ‘on this observation, we defined a standing molecule as one for which the C=C axis tilt exceeds 65’) with respect to the sur-
35
55
75
95
T (K> Fig. 4. The temperature variation of (a) the molecule-molecule (m-m) and (b) the molecule-surface (m-s) contributions to the average configurational energy. E.
face plane Fig. 2b shows a plot of the fraction of such molecules versus temperature. This quantity increased steadily in the liquid phase. In the solid phase prior to melting, random molecules were observed to stand erect and stay that way for times of order lo20 ps. We also always observed dislocations in the solid phase prior to melting. However, in view of the small size of the simulation systems we employed, it was not possible to establish whether or not fluctuations involving the standing up of molecules were unequivocally correlated with the formation and motion of dislocations which, in turn, drive the melting, as in the Kosterlitz-Thouless theory [ 111. From an inspection study of the phase space trajectories, it was evident that local density fluctuations involving standing molecules indicated the onset of melting in the system with n=0.84 (see fig. 2b). In the other system (n=0.90), these fluctuations involving standing molecules persist for longer times and on average, there are already a few molecules standing in the rotator phase solid. These we attribute to the fact that the II = 0.90 system is in a compressed state, since at high temperature ( T> 40 K), our model has
M.A. Moller, M.L. Klein /The melting of ethylene on graphite
an equilibrium density of 0.84 [ 41. In this case, the fluctuations serve to relieve the stress, but the system remains solid. From this we conclude that the standing of molecules and the formation of associated dislocations is not alone sufficient to drive the melting. However, the melting process seems to be the same in the low-density case (cf. heat capacity anomaly in fig. 4b). Fig. ,l b shows the evolution of the structure factor as the n = 0.90 system goes through the melting transition. The number of standing molecules was monitored throughout this run. The fraction of erect molecules varies smoothly from 0.05 in the solid region ( t < 100 ps) to 0.10 in the liquid region ( t > 100 ps). The steady change in the number of erect molecules as the system proceeds through the melting transition likely does not drive the transition, but acts to smear it, and hence make it appear to be continuous. In any case, fig. 4 confirms that our simulated transition definitely has a broad thermal anomaly at both coverages, in excellent agreement with experiment. No second layer promotion of molecules is observed in the melting process.
239
at the experimentally determined surface density would be expected to show enhanced rotational and translational motion in the solid. Neither of these points will likely have any bearing at all on the nature of the melting transition, which we conclude is intimately related to the gradual standing up of the ethylene molecules. A similar phenomenon has recently been observed in MD calculations on the system O,/ graphite [ 13 1.
Acknowledgement We thank Moses Chair, Toni Heidemann, John Larese, Larry Passell, Dieter Richter, Bill Steele, and especially the late Jim Morrison for much helpful discussion and correspondence. This work was supported in part by the National Sciences and Engineering Research Council of Canada.
References [ I ] H.K. Kim and M.H.W. Chan, Phys. Rev. Letters 53 (1984)
4. Discussion We have presented the result of a molecular dynamics simulation of the melting of ethylene on a corrugated graphite surface. The melting temperature of the computer model is lower than the experimental value. Two possible reasons for this come to mind. Firstly, we have not taken into account the anisotropy of the molecule-surface potential, which is due to the greater polarizability of the graphite within the plane of the surface rather than normal to it. This anisotropy has been shown to give rise to an increased barrier to translation across the plane of the surface. In the case of NZ, this lead to an increase of 9” in the melting temperature for a similar computer model [ 12 1. Secondly, there are likely residual errors in the intermolecular potential itself. In earlier work, this potential gave solid structures which were too dense [ 41. Due to this fact, a simulation carried out
170. [2]H.K. Kim, Q.M. Zhang and M.H.W. Chat-r, Phys. Rev. Letters 56 (1986) 1579; A.D. Migone, Z.R. Li and M.H.W. Char-r, Phys. Rev. Letters 53 (1984) 810. [ 31 E.D. Specht, R.J. Birgeneau, K.L. D’Amico, D.E. Moncton, S.E. Nagler and P.M. Horn, J. Phys. Lett. 46 (1985) L561; SE. Nagler, P.M. Horn, T.F. Rosenbaum, R.J. Birgeneau, M. Sutton, S.G.J. Mochrie, D.E. Moncton and R. Clarke, Phys. Rev. B 32 (1985) 7373. [4] S. Nose and M.L. Klein, Phys. Rev. Letters 53 ( 1984) 818. [5] M.L. Klein and J.A. Morrison, Carbon 25 (1987) 23. [6]‘D.E. Williams, J. Chem. Phys. 45 (1966) 3770. [ 71 W.A. Steele, Surface Sci. 36 ( 1973) 317. [ 81 B.H. Grier, L. Passell, J. Eckert, H. Patterson, D. Richter and R.J. Rollefson, Phys. Rev. Letters 53 ( 1984) 8 14; S.K. Satija, L. Passell, J. Eckert, W. Ellenson and H. Patterson, Phys. Rev. Letters 5 1 ( 1983 ) 4 11. [9] J.Z. LareseandR.J. Rollefson, Phys. Rev. B 31 (1985) 3048. [ 10 ] J.Z. Larese, L. Passell, A.D. Heidemann, D. Richter and J.P. Wicksted, Phys. Rev. Letters 61 ( 1988) 432. [ 1 1] J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6 ( 1973) 1181. [ 121 Y.P. Joshi and D.J. Tildesley, Mol. Phys. 55 ( 1985) 999. [ 131 V.R. Bhethanabotla and W.A. Steele, Langmuir 3 ( 1987) 581.
235-239
THE CONTINUOUS MELTING TRANSITION A MOLECULAR DYNAMICS STUDY
OF ETHYLENE
ON GRAPHITE:
Michael A. MOLLER Department of Chemistry, McMaster University, Hamilton, ‘Ontario. Canada LBS 4Ml
and Michael L. KLEIN of Chemistry,University ofPennsylvania, Philadelphia, PA 19104-6323, USA
Department Received
1 May 1988
The melting of monolayer films of ethylene on graphite has been investigated by molecular dynamics calculations. Our results suggest that melting is influenced by the presence of individual molecules standing erect, with their C=C axes perpendicular to the graphite surface. These molecules lower the effective density in their vicinity and because their number increases with temperature, cause the melting to appear continuous.
1. Introduction
graphite identified translation-rotation coupling as playing an important role in the dynamics of ethylene overlayers [ 41. However, the process in which molecules become erect does not appear to be related to such cooperative (collective) phenomena or a phase separation. The fluctuations seem to occur randomly in the sample, and persist for times of order lo-20 ps. We suggest that the observed anomalous melting transition for ehtylene on graphite is intimately associated with this novel process. We now describe the MD calculations which lead us to this conclusion.
The diversity of the melting characteristics exhibited by quasi two-dimensional molecular overlayers, physisorbed on graphite, is quite surprising [ l-3 1. For example, submonolayer films of the simplest hydrocarbon, methane (CH,), show a sharp first-order transition; the heat capacity peak associated with melting is narrow and exceeds 100 k, in height [ 11. The behavior of the aspherical molecule ehtylene (C,H,) is quite different. Here, melting extends over several degrees and the heat capacity peak has a maximum height of only about 2 k, [ 21. The latter observation has been interpreted as evidence for a continuous melting transition [ 2 1. In order to seek a possible microscopic explanation of this observation, we therefore embarked on a molecular dynamics (MD) study of the ethylene-graphite system. An analysis of our MD calculations leads us to conclude that the melting of ethylene films physisorbed on graphite is strongly influenced by molecular reorientation. The local density fluctuations required for melting are enhanced when some of the molecules stand up, with their C=C axes perpendicular to the graphite basal plane. A previous MD study of C,H,/ 0301-0104/89/$ ( North-Holland
03.50 0 Elsevier Science Publishers Physics Publishing Division )
2. Calculations We have employed the same model as used previously to discuss the structure and dynamics of the solid phases formed by ethylene on graphite at low temperatures [ 4,5 1. The intermolecular potential for the ethylene molecules was the sum of atom-atom interactions of the exp-6 type (see ref. [ 61, “the potential set IV”). The molecule-surface potential was derived from atom-atom interactions between the atoms of the ethylene molecules and the atoms comB.V.
236
M.A. Moller. M.L. Klein /The melting ofethylene on graphite
prising the graphite basal plane. In order to include the corrugation, we have followed the method of Steele [ 7 1. Here, .the atom-surface potentials are developed in a Fourier series expansion. The first term, V,,(z), yields the flat surface approximation used in our earlier MD studies of this system [4]. The second term in the Fourier series, V, (x,y,z), contains the effect of the surface corrugation. The variation of V,, as an ethylene molecule moves across the graphite surface is about 50 K for the potentials used here. The optimal configurational energy for a single molecule interacting with the graphite basal plane agrees to within a few percent of the value derived from the measured isosteric heat of adsorption extrapolated to zero coverage, namely - 19.9 kJ/mol [ 51. Inclusion of the corrugation term V, imposes an important constraint on the system to be studied. Since the simulation box had periodic boundaries parallel to the (x,y) plane, the potential had to be continuous at these boundaries. To ensure that this was so, the boundaries were kept fixed, with an integral number of graphite hexagons inside the MD box. Simulations were performed in the (N,E, V) ensemble at two different coverages. We define a coverage of y1= 1 to be a commensurate &X v1? R30” monolayer. One set of calculations employed a box with sides 13$a x 23a, where a= 2.4 A is the lattice constant for the graphite basal plane hexagons. There were 598 hexagons within this box. We used N= 168 ethylene molecules, so the coverage in this case was n = 0.84. This value corresponds to the coverage of a complete monolayer of the so-called low-density phase at low temperature. The remaining simulations were carried out with N= 132 molecules in a box with edges 11 ,/?a x 20a. Since this box contained 440 hexagons, the coverage was n = 0.90. This was the lowtemperature equilibrium coverage for our potential model as determined in previous calculations [ 41. By running simulations at this coverage, we eliminate the possibility that the low melting temperature at n=0.84 and the resulting enhanced motion of the molecules is due to a coverage less than the equilibrium value for our model. Experimentally, at low temperature and at coverages less than 0.84, solid ethylene overlayers adopt an incommensurate structure, expanded with respect to the registered phase, in which the C=C molecular axes are predominantly parallel to the graphite basal plane
[4,8]. Moreover, before it melts, this so-called lowdensity solid transforms (at Tz 35 K) to a rotator phase [ 2,4,8,9]. In view of these facts, our MD calculations were initiated from ordered triangular herringbone structures with all C=C axes parallel to the graphite surface. Since we are interested in melting, systems were prepared at temperatures ranging from 40 to 100 K [ 8 1. In all cases, the ethylene molecules became orientationally disordered in the plane of the surface. A typical simulation, run for about 200 ps using a time step of 2.5 fs, took about 20 h of CPU time on an IBM 3090; it was therefore impractical to investigate a larger system size. To locate the melting temperature, the structure amplitude
S(k)+
cos(k.r))
was evaluated as a running average throughout each simulation. Here, k is one of the three smallest reciprocal lattice vectors of the triangular lattice appropriate to the surface coverage of interest, Yis the distance between the centers of mass of two molecules and the summation is over all such pairs of molecules. The structure amplitude was calculated for each molecule in the system as the origin, and then averaged. The average was also taken over all pairs of molecules and time steps. As defined, S(k) = 1 for a perfectly static triangular lattice and falls to near zero for a translationally disordered liquid. In a solid at temperatures above zero, S(k) is reduced to a value less than unity by the Debye-Waller effect (thermal motions), and static disorder (dislocations). In the following discussion and figures, we have taken the structure amplitude to be the average of all three S(k) values where k=(l,O), (0, l)and(l, 1). We also studied in detail the diffusional motions of the ethylene molecules on the surface. The translational diffusion coefficient, D,, was obtained from the time dependence of the mean square displacements of the molecular centers of mass in the (x, Ji) plane. Rotational diffusion coefficients, D,,, were obtained from the appropriate angular velocity autocorrelation functions. We followed rotation about axes in the molecular plane both parallel and perpendicular to the C=C axis as well as rotation about an axis perpendicular to the molecular plane. Rota-
M.A. Moller, M.L. Klein /The melting of ethylene on graphite
tional motion about C=C is by far the most important. Other quantities monitored included the average tilt of the C=C axis away from the plane of the surface, and the configurational energy, which was separated into a molecule-molecule part E( m-m) and a molecule-surface part E( m-s).
0.6 0.5
ti! 3
0.4
6
0.3
E =
0.2
I I
_I
3. Results
0.0
%
0.6
5
O*
E
0.2 0.0 0.6
2
04 .
5 0.2
0
50
toot5o200250
TM @> Fig. 1. Time evolution of the structure amplitude for ethylene overlayers physisorbed on graphite. The upper curve (a) is for T= 53 K and the lower curve (b) is for T= 63 K. Both calculations were carried out at a surface coverage of n=0.90 (see text).
I
l
I
vO.64
. n=O.90
t
,I
*fl t,,
I
’
b
I
, ’
, * I .
- 04 0.15
t
t:
t
is ii?
8
II
t7; 0.x)
&
I
t (4
9 B 3
8
I’t
0.1 _
Fig. 1 shows the time evolution of the structure amplitude S(k) for points on either side of the melting transition in the system with n=0.90. Close to the melting transition, long simulation trajectories (200 ps) are required to establish equilibrium. In fig. 2a we show the temperature dependence of S( k) for both coverages. The value of S(k) was calculated every 25 steps, and these values were averaged to produce the figure. The error bars represent the standard deviation of S(k). For the series with n=0.84 (dots) there is a decrease in S(k) in the region 5260 K whereas the drop occurs between 57 and 63 K for the series with n = 0.90 (squares). This drop-off
231
t
0.05
t tt + +*
tt I
0.00 35
55
75
95
TEw09 Fig. 2. (a) The temperature dependence of the structure amplitude and (b) the fraction of molecules in the system standing up.
in S(k) identified the melting region for a given coverage. The fluctuations increased considerably as the temperature approached the melting point, representing the formation and annealing of dislocations. Unfortunately, the calculated melting regions are somewhat lower than the experimental counterparts which are centered around 70 and 80 K. This disagreement likely signals inadequacies in the relevant interaction potentials we have employed. Fig. 3 shows the temperature variation of the diffusion coefficients D,, and D,,. A line drawn through the values of D*, extrapolates to zero at about 50 K for n~0.84 and at about 55 Kfor n~0.90. These values lie below the respective melting regions and suggest that D,, is continuous across the melting “transition”. Recent quasi-elastic neutron scattering data are in agreement with this result but also suggest that the melting is associated with enhanced molecular rotation about the molecular axes perpendicular to C=C [ lo]. This latter observation seems to be at variance with our MD results (fig. 3b). This discrepancy may be related to the fact that, in the experiment, the overlayer can expand on heating, whereas the calculations are carried out at fixed coverage, which would restrict the motions of the molecules.
238
M.A. Moller, M.L. Klein /The melting qfethplene
c
I
’
r
I
I
on graphite
-6.2
(a) . n=O.84
w,
t
. .
(al
. n=O.90
-6.9
-7.6
-y rn
9.0
Fg
6.0
0
P
t
t
+k
t
tt
t+:
3.0
+ t t I
0.0 35
A
Ati
b
I
55
T
(;
@*)! @*I 4
I
95
Fig. 3. (a) Temperature variation of the translational diffusion coefficient, D,,. (b)Temperature variation of the rotational diffusion coefficients, D,,,, for nE0.84. Symbols indicate rotation about the C-C axis (circles), and about orthogonal axes perpendicular (triangles), and parallel (crosses) to the molecular plane.
Fig. 4 shows the temperature variation of the configurational energies E(m-m) and E(m-s). These quantities exhibit only a small variation across the melting region. Indeed, crude estimates of the derivatives of these curves with respect to temperature suggest a specific heat anomaly of order 2-3 k, in the relevant transition regions. This behavior is entirely consistent with the measured heat capacity data [ 21. Even though the calculated melting temperature is somewhat low, there is apparent success in modeling the thermal anomaly at the melting transition [ 21. Accordingly, we were encouraged to search for a possible microscopic explanation of this weak transition. Of all the quantities we examined, only one showed any dramatic change at melting and that was the number of molecules standing up in the system. The orientational distribution function for C=C with respect to the surface plane was bimodal, with peaks at 0” and 80” and a well-defined minimum extending from about 20” to 65”. Based ‘on this observation, we defined a standing molecule as one for which the C=C axis tilt exceeds 65’) with respect to the sur-
35
55
75
95
T (K> Fig. 4. The temperature variation of (a) the molecule-molecule (m-m) and (b) the molecule-surface (m-s) contributions to the average configurational energy. E.
face plane Fig. 2b shows a plot of the fraction of such molecules versus temperature. This quantity increased steadily in the liquid phase. In the solid phase prior to melting, random molecules were observed to stand erect and stay that way for times of order lo20 ps. We also always observed dislocations in the solid phase prior to melting. However, in view of the small size of the simulation systems we employed, it was not possible to establish whether or not fluctuations involving the standing up of molecules were unequivocally correlated with the formation and motion of dislocations which, in turn, drive the melting, as in the Kosterlitz-Thouless theory [ 111. From an inspection study of the phase space trajectories, it was evident that local density fluctuations involving standing molecules indicated the onset of melting in the system with n=0.84 (see fig. 2b). In the other system (n=0.90), these fluctuations involving standing molecules persist for longer times and on average, there are already a few molecules standing in the rotator phase solid. These we attribute to the fact that the II = 0.90 system is in a compressed state, since at high temperature ( T> 40 K), our model has
M.A. Moller, M.L. Klein /The melting of ethylene on graphite
an equilibrium density of 0.84 [ 41. In this case, the fluctuations serve to relieve the stress, but the system remains solid. From this we conclude that the standing of molecules and the formation of associated dislocations is not alone sufficient to drive the melting. However, the melting process seems to be the same in the low-density case (cf. heat capacity anomaly in fig. 4b). Fig. ,l b shows the evolution of the structure factor as the n = 0.90 system goes through the melting transition. The number of standing molecules was monitored throughout this run. The fraction of erect molecules varies smoothly from 0.05 in the solid region ( t < 100 ps) to 0.10 in the liquid region ( t > 100 ps). The steady change in the number of erect molecules as the system proceeds through the melting transition likely does not drive the transition, but acts to smear it, and hence make it appear to be continuous. In any case, fig. 4 confirms that our simulated transition definitely has a broad thermal anomaly at both coverages, in excellent agreement with experiment. No second layer promotion of molecules is observed in the melting process.
239
at the experimentally determined surface density would be expected to show enhanced rotational and translational motion in the solid. Neither of these points will likely have any bearing at all on the nature of the melting transition, which we conclude is intimately related to the gradual standing up of the ethylene molecules. A similar phenomenon has recently been observed in MD calculations on the system O,/ graphite [ 13 1.
Acknowledgement We thank Moses Chair, Toni Heidemann, John Larese, Larry Passell, Dieter Richter, Bill Steele, and especially the late Jim Morrison for much helpful discussion and correspondence. This work was supported in part by the National Sciences and Engineering Research Council of Canada.
References [ I ] H.K. Kim and M.H.W. Chan, Phys. Rev. Letters 53 (1984)
4. Discussion We have presented the result of a molecular dynamics simulation of the melting of ethylene on a corrugated graphite surface. The melting temperature of the computer model is lower than the experimental value. Two possible reasons for this come to mind. Firstly, we have not taken into account the anisotropy of the molecule-surface potential, which is due to the greater polarizability of the graphite within the plane of the surface rather than normal to it. This anisotropy has been shown to give rise to an increased barrier to translation across the plane of the surface. In the case of NZ, this lead to an increase of 9” in the melting temperature for a similar computer model [ 12 1. Secondly, there are likely residual errors in the intermolecular potential itself. In earlier work, this potential gave solid structures which were too dense [ 41. Due to this fact, a simulation carried out
170. [2]H.K. Kim, Q.M. Zhang and M.H.W. Chat-r, Phys. Rev. Letters 56 (1986) 1579; A.D. Migone, Z.R. Li and M.H.W. Char-r, Phys. Rev. Letters 53 (1984) 810. [ 31 E.D. Specht, R.J. Birgeneau, K.L. D’Amico, D.E. Moncton, S.E. Nagler and P.M. Horn, J. Phys. Lett. 46 (1985) L561; SE. Nagler, P.M. Horn, T.F. Rosenbaum, R.J. Birgeneau, M. Sutton, S.G.J. Mochrie, D.E. Moncton and R. Clarke, Phys. Rev. B 32 (1985) 7373. [4] S. Nose and M.L. Klein, Phys. Rev. Letters 53 ( 1984) 818. [5] M.L. Klein and J.A. Morrison, Carbon 25 (1987) 23. [6]‘D.E. Williams, J. Chem. Phys. 45 (1966) 3770. [ 71 W.A. Steele, Surface Sci. 36 ( 1973) 317. [ 81 B.H. Grier, L. Passell, J. Eckert, H. Patterson, D. Richter and R.J. Rollefson, Phys. Rev. Letters 53 ( 1984) 8 14; S.K. Satija, L. Passell, J. Eckert, W. Ellenson and H. Patterson, Phys. Rev. Letters 5 1 ( 1983 ) 4 11. [9] J.Z. LareseandR.J. Rollefson, Phys. Rev. B 31 (1985) 3048. [ 10 ] J.Z. Larese, L. Passell, A.D. Heidemann, D. Richter and J.P. Wicksted, Phys. Rev. Letters 61 ( 1988) 432. [ 1 1] J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6 ( 1973) 1181. [ 121 Y.P. Joshi and D.J. Tildesley, Mol. Phys. 55 ( 1985) 999. [ 131 V.R. Bhethanabotla and W.A. Steele, Langmuir 3 ( 1987) 581.