Langmuir 1987,3, 581-587
581
Molecular Dynamics Simulations of Oxygen Monolayers on Graphite Venkat R. Bhethanabotla and William A. Steele* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802 Received September 23, 1986. In Final Form: February 13,1987 Molecular dynamics computer simulations are reported for oxygen monolayers adsorbed on the graphite basal plane. Intermolecular potentials were modeled by using a site-site description that represents both the nonspherical shape for the O2and the periodic variations in the 02-solid interactions. Systems studied included a complete layer with an initial configuration taken to be close to that of the experimental &phase and a strip of O2 molecules initially covering half the surface, with a &phase configuration in the strip. Temperatures covered the range from the solid phase at -20 K through the melting transition located at 27 and 36 K for the patch and the full surface, respectively, and continuing up to -70 K, which is above the two-dimensional liquid-vapor critical point for the patch system. Translational and orientational order parameters were determined, together with time-dependent trajectories and associated functions which allow one to characterize the dynamical properties of the adsorbed molecules. The general picture that emerges is a simultaneous loss of translational and orientational order at the melting points. Sliding of the solid layer across the adsorbent surface was observed for the full surface system. The sudden alteration in the out-of-plane orientational properties observed at the melting point for the full surface system is believed to play a significant role in the melting process for this phase.
Introduction I t is now known that monolayers physisorbed on the uniform graphite basal plane form a variety of ordered quasi-twedimensional solid phases. When the physisorbed molecules are nonspherical, the phases are orientationally as well as translationally ordered a t low temperatures. As a consequence, the phase diagrams of such systems can become quite complex. Possibly the most extensively investigated system of this kind is nitrogen. A combination of diffraction and thermodynamic measurements1s2indicated that the commensurate phase undergoes a transition at -27 K from herringbone orientationally ordered to an orientationally disordered phase. This phase then melts directly to a supercritical fluid a t 45 K with no intervening two-dimensional liquid. At higher surface densities, a t least two incommensurate solids are formed which are orientationally ordered a t low temperature in slightly altered herringbone structures relative to that for the commensurate phase. For the most part, these experimental findings have been confirmed and extended by computer
sir nu la ti on^.^^^ (1)Kjem, J. K.; Passell, L.; Taub, H.; Dash, J. G.; Novaco, A. D. Phys. Rev. B 1976,13,1446.Eckert, J.;Elleneon, W. D.; Hastings, J. B.; Passell, L. Phys. Rev. Lett. 1979,43,1329. Diehl, R.D.; Toney, M. F.; Fain, S. C., Jr. Phys. Reu. Lett. 1982,423,177. Diehl, R.D.; Fain, S. C., Jr. Surf. Sci. 1983,125,116.Fain, S.C., Jr.; Chinn, M. D.; Diehl, R. D. Phys. Reu. B 1980,21,4170. Nielsen, M.; Kjaer, K.; Bohr J . Electron Spectrosc. Relat. Phenom. 1983,30,111. Morishige, K.; Mowforth, C.; Thomas, R. K. Surf. Sci. 1985,151,289.Diehl, R.D.; Fain, S. Phys. Rev. B 1982,26, 4785. Diehl, R.D.; Fain, S. C., Jr. J. Chem. Phys. 1982,77,5065. You, H. Fain, S. C., Jr. Discuss. Faraday SOC.,in press. Fain, S. C. Ber. Bumenges. Phys. Chem. 1986,90,211. (2)Butler, D. M.; Huff, G. B.; Toth, R. W.; Stewart, G. A. Phys. Rev. Lett. 1975,35,1718.Migone, A. D.;Kim, H. K.; Chan, M. H. W.; Talbot, J.; Tildesley, D. J.; Steele, W. A. Phys. Rev. Lett. 1983,51, 192. Zhang, Q.M.; Kim, H. K.; Chan, M. H. W. Phys. Rev. Lett. B 1986,33E,413. Zhang, Q. M.; Kim, H. K.; Chan, M. H. W. Phys. Rev. B 1985,32,1820. Chan, M. H.W.; Migone, A. D.; Miner, K. D.; Li, 2.R. Phys. Rev. B 1984, 30,2681. Piper, J.; Morrison, J. A.; Peters, C.; Ozaki, Y. J. Chem. SOC., Faraday Trans. 1983,79,2533.Chung, T. T.; Dash, J. G. Surf. Sci. 1977, 66, 559. (3)Talbot, J.; Tildesley, D. J.; Steele, W. A. Mol. Phys. 1984,51,1331. Vernov, A. V.; Steele, W. A. Langmuir 1986,2,219.Vernov, A. V.; Steele, W. A. Surf. Sci. 1986,171,83. Peters, C.; Klein, M. Mol. Phys. 1985,54, 985. Joshi, Y. P.; Tildesley, D. J. Mol. Phys. 1985,55,999. Talbot, J.; Tildesley, D. J.; Steele, W. A. Discuss. Faraday SOC.,in press. Vernov, A. V.; Steele, W. A. Langmuir 1986,2,219.
0743-7463I87 12403-0581$01.50,IO ,
I
Extensive experimental studies of the properties of physisorbed oxygen monolayers on graphite have also been reported."' In spite of the general similarity in size and shape of O2and N2, the phase diagrams of the two adsorbed systems are quite different.5 In part, this is due to the presence of magnetically ordered O2phases at very low temperatures. However, even at temperatures ranging from -15 to -55 K where the magnetic part of the interaction no longer plays a significant role, ordered O2 monolayers are distinctly different from those of N2. In the first place, the commensurate translational ordering which is so often observed in simple solid monolayers on graphite (N2, Kr, Xe, CHI, etc.) is absent in 02.Furthermore, the herringbone orientational ordering observed in N2,for which 4/6 of the nearest-neighbor pairs are nearly perpendicular to one another, is absent in the O2 case. In fact, the low-density ordered phase in O2 is an incommensurate array of centered rectangular cells with orientational ordering in which all molecular axes are parallel to one another (6 phase). This phase can be compressed, but the low-densityform melts directly to a liquid a t -26 K without showing an analogue to the orientationally disordered solid observed in N2 layers a t T > 29 K. Adsorption isotherms indicate liquid-gas coexistence with a two-dimensional critical temperature of -65 K. For convenience, a portion of the experimental phase diagram (4)Talbot, J.; Tildesley, D. J.; Steele, W. A. Surf. Sci. 1986,169,71. (5)Fain, S.C., Jr.; Toney, M. F.; Diehl, R. D. In Proceedings of Ninth International Vacuum Congress and Fifth International Conference on Solid Surfaces;desegovia, J. L. Ed.; Imprenta Moderna: Madrid, 1983; p 129. Toney, M. F.;Fain, S. C., Jr., submitted for publication in Phys. Reu. E . (6)Toney, M. F.; Diehl, R. D.; Fain, S.C., Jr. Phys. Rev. B 1983,27, 6413. McTague, J. P.; Nielsen, M. Phys. Reu. Lett. 1976,37,596.Nielsen, M.; McTague, J. P. Phys. Rev. B 1976,19,3096.Heiney, P. A.; Stephens, P. W.; Mochrie, S. G.; Akimitau, J.; Birgeneau, R. Surf. Sci. 1983, 125, 539. Stephens, P. W.;Heiney, P. A.; Birgeneau, J. R.; Horn, P. M.; Stoltenberg, J.; Vilches, 0. E. Phys. Rev. Lett. 1980,45,1959. Mochrie, S. G.; Sutton, M.; Akimitau, J.; Birgeneau, R. J.; Horn, P. M.; Gregory, S. Phys. Rev. Lett. 1978,40,723;1977,39,1035. Awschalom, D.; Lewis, G.; Gregory, S. Phys. Rev. Lett. 1983,51,586. Stoltenberg, J.; Vilches, 0. Phys. Rev. B 1980,22,2920.Marx, R.;Christoffer, B. Phys. Reu. Lett. 1983,51,790.Marx, R.;Braun, R. Solid State Comrnun. 1980,33,229. Yon, H.; Fain, S. C., Jr. Phys. Rev. B 1986,33,5886. (7)Dericbourg, J. Surf. Sci. 1976,59,554. Gilquin, J. Ph.D. Thesis, Saclay, 1979.
0 - 1987 American Chemical Societv
582 Langmuir, Vol. 3, No. 4 , 1987 8.5.
I
I
; BtS, -
I
I
i
I‘ + 3 ,
Bhethanabotla and Steele
+( I
1
I
Table I. Site-Site Potential Parameters” e l k , K u, A tlk, K 0-0 52.1 2.99 N-N 36.3 0-carbon 38.2 3.195 N-carbon 31.9 28.0 3.40 carbon-carbon
t fluid
I
f
u, A
3.31 3.36
Bond length = 1.21 A for 02;1.10 A for N,.
IO
I
I
15
20
II
25
I
I
30 35 T (K)
I
1
40
45
1
50
Figure 1. Portion of the experimentally determined phase diagram5for O2 on gcb is reproduced here. Circles show the temperatures and coverages of the simulations reported here.
for O2 on graphite is reproduced in Figure 1. In view of the significant differences between O2 and N2 monolayer phase behavior, it is obviously of interest to perform comparative calculations on the two systems. A number of such studies have already been reported for 02: energy minimi~ations,8,~ which yield structures in the limit of 0 K:, lattice dynamics calculations of the vibrational and librational frequencies in O2 1ayers;’O and calculations of magnetic 0rdering.l’ In this paper, molecular dynamics simulations of an O2 monolayer adsorbed on the graphite basal plane are reported. Intermolecular potential functions are used that give a good representation of the known second bulk virial coefficient and bulk liquid and solidstate properties.12 In this way, the shape, size, well-depth, and the electrostatic part of the interactions are constrained to fit other experimental data. The 02-solid potential was constructed by analogy with the N2 and the rare-gas molecule-solid interactions used successfully in other calculations.3J3J4 Simulations were performed at temperatures ranging from 10 to 75 K; initial densities (and monolayer structure) were taken to be close to the known low-temperature G.phase values. Two different cases were studied: one in which the monolayer filled the surface, so that density within the layer was constrained to be constant at all temperatures (in the absence of promotion of a molecule to the second layer), and another set of simulations of a “patch” in which the initial configuration of molecules covered only half the surface and thus could expand (or evaporate) to a lower density equilibrium phase. Phase points studied in these simulations are indicated in Figure 1. We find that the intermolecular interaction laws used here are capable of reproducing the experimental data rather well. In addition, the detailed molecular information generated in the simulations allows us to give a new and more detailed interpretation of the molecular processes involved in the phase transitions for these systems. (8) Etters, R. D.; Pan, R. P.; Chandrasekharan, V. Phys. Reu. Lett.
1980,45,645. Pan, R. P.; Etters, R. D.; Kobashi, K.; Chandrasenkharan, V. J. Chem. Phys. 1982, 77, 1035. Flurchick, K.; Etters, R. D. J . Chem. Phys. 1986,84, 4657. (9) Joshi, Y. P.; Tildesley, D. J. Surf. Sci. 1986, 166, 169. (10) Etters, R. D.; Kobashi, K. J. Chem. Phys. 1984, 81, 6249. (11) Duparc, 0. M. B.; Etters, R. D. J. Chem. Phys. 1986, 86, 1020. (12) Klein, M. L.; Levesque, D.; Weis, J. J. Phys. Rev. 1980,21, 5785. (13) Steele, W. A. J . Phys. Chem. 1978, 82, 817. (14) Bojan, M. J.; Steele, W. A. Langmuir 1987, 3. 116.
Simulations The simulation technique used was not significantly different from that employed in previous studies of adsorbed N2.3 Molecular motions were calculated with a conventional molecular dynamics15 predictor-corrector algorithm16using periodic boundary conditions in the plane of the surface. Some trajectories were evaluated in a constant energy (microcanonical) ensemble. Since a molecule could occasionally escape from the surface at high temperature, a constant-temperature algorithm1’ is more appropriate, since it avoids the large temperature fluctuations produced in the microcanonical ensemble when a molecule undergoes a large potential energy change. In fact, most of the simulations reported here are for constant N,V,T.l8 A finite time step of length 5 X lov3ps was used throughout in the integration of the equations of motion. (This gave excellent energy conservation in the N ,V,E simulations.) As it turned out, promotion to the second layer was not observed except for T 2 60 K. Quaternions were used to represent the orientations of the rigid diatomic molequles and a number of time-dependent translational and rotational properties were evaluated as well as the relevant thermodynamic data such as the component parts of the average potential energy of an adsorbed molecule. Samples of 192 O2 molecules on a surface of appropriate area were simulated; in the “patch” studies, the initial configuration was taken to extend from boundary to boundary in one direction, so that the periodic boundary conditions generated infinitely long strips of adsorbed solid separated by a vacant surface into which the molecules could “evaporate” or diffuse, if the temperature was sufficiently high. The simulations were first run for intervals of 3000-4000 time steps to achieve equilibrium. Data were gathered subsequently for intervals ranging upward from 2000 time steps; positions and velocities at intervals of 10 time steps were saved for future analysis, but thermodynamic averages were calculated during the run. Table I shows the parameters for the intermolecular potentials that were used in this work. The 02-02 interaction is represented by a two-site Lennard-Jones model in which spherically symmetric interactions were assumed between two sites located at the 0 atom positions within a molecule and the pairs of similar sites on neighboring O2 molecules. Well depths and size parameters for the 12-6 sitesite function are listed in Table I-these were adjusted to give agreement between calculated bulk properties such as solid-state structural and dynamic properties or second virial coefficient and experiment. (A limited number of bulk liquid simulations were also carried out by us-these also compared favorably with experiment.) Note that both the magnetic and the electrostatic quadrupolar interactions (15) See, for example: Streett, W. B.; Tildesley, D. J. In Computer Modelling of Matter; Lykos, P., Ed.; ACS Symposium Series 86; Murad, S.; Gubbins, K. In Computer Modelling of Matter; Lykos, P., Ed.; ACS Symposium Series 86; American Chemical Society: Washington, DC, 1976. Evans, D. J. Mol. Phys. 1977, 34, 317. (16) Gear, C. W. In Numerical Initial Value Problems in Ordinary Differential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1971. (17) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, I , 297; Chem. Phys. 1984,87, 451. (18) Nos& S. J . Chem. Phys. 1984, 81, 511.
Langmuir, Vol. 3, No. 4 , 1987 583
Simulations of Oxygen Monolayers on Graphite
are quite small for O2 and were omitted from the simulation on the grounds that their effects on the properties of these systems would be negligible at the temperatures studied. The 02-graphite interaction was modeled by using an extension of the site-site representation in which the graphite is taken to be an ordered array of C sites, each of which interacts with an O2 via spherically symmetric 0 s i t e 4 site Lennard-Jones functions. Parameters of this function were obtained by applying Lorentz-Berthelot combining rules (geometric mean for well depth c/k and arithmetic mean for size a). The C-C “site” parameters were taken from previous work on the N2/graphite and rare gas/graphite systems. The infinite sum over C sites needed to obtain the 0 site-solid interaction potential was expressed as a Fourier series with the periodicity of graphite basal plane and the nonnegligible coefficients were evaluated by using the prescription of Steele.lg The most straightforward test of the model moleculesolid potential is to compare the simulated value of ( cgs), the average molecule-solid potential, with the experimental heat of adsorption qst, extrapolated to zero coverage, minus RT. Although there is a minor problem of matching the experimental and the simulation temperature, the data give qst(0)/R- T = 1156 K6J4 (with an uncertainty of several percent) which agrees with the simulation value of 1164 K, to well within the combined uncertainties of the two numbers. One interesting aspect of these simulations is the comparison and contrast between the properties of adsorbed O2 and N2. Thus, the parameters for the site-site model for NZ3J4are listed in Table I together with those for 02. The first point to note is that the Lorentz-Berthelot combining rules applied to known site-site potentials for the bulk fluid and the carbon-carbon parameters listed in Table I yield molecular-solid potentials that reproduce both the experimental qJ0) and the low-coverage limiting free energies (measured via the Henry’s Law constants)’* to quite high accuracy. This is strong evidence of the utility of this approach for molecules on graphitized carbon black (gcb) (and similar solids). The rare gas/gcb potentials appear to be accurately parameterized in this way,13 but the extension of the basic idea to molecular adsorbates would be a convenient and widely applicable technique, if it can be shown to be accurate. Second, it should be noted that the energy minimization calculations give quantitative support to the observation that the O2 molecule is somewhat smaller than N2, with approximate width and length parameters of 3.0 and 4.2 8, compared to 3.3 and 4.4 8, for N2. Consequently, the formation of an ordered phase for adsorbed O2with density greater than the commensurate value observed for N2 is expected. Furthermore, it is the electrostatic quadrupole-quadrupole energy which dictates the herringbone pattern observed for N2. In its absence, it is reasonable to find O2in a structure with all molecular axes parallel. The remaining question which we hope to answer in this paper is: what is the effect of increasing temperature on this two-dimensional O2crystal?
Equilibrium Properties In a computer simulation, it is convenient to evaluate the average energy per molecule a t each temperature considered. Of course, the average kinetic energy for a (19)Steele, W. A. Surf. Sci. 1973, 36, 317; J.Phys. (Paris) 1977, C4, 61. (20) Sacco, J.
5071.
E.;Sokoloff, J. B.;Widom, A. Phys. Rev. B
1979, 20,
Patch 0
Full surface
-1
- 3 4
$4 p
-31
1 Temperature
Figure 2. Simulated values of the average molecule-solid potential energy (in reduced units) are plotted vs. temperature. Reduced units are energylc,,, where t,,lk = 38.2 K.
\I
-7-5i
Ib
o:
io
20
510
do
Ternperat ure
Figure 3. Average 02-02interaction energy is plotted vs. temperature for the two systems considered. T h e full surface is essentially a constant-density system, but the density in the patch decreases with thermal expansion, with melting a t 27 K, and, eventually, with evaporation. No step for evaporation is observed for T between 60 and 70 K. This is not unexpected, since the calculations are sparse in this region and the density is close to critical, which minimizes AE for vaporization.
linear molecule is just 5/2kT,so that the quantity of interest is the average potential energy. For adsorbed monolayers, one can evaluate the components of this energy which are, for O2 on gcb, (egs), the average molecule-solid interaction mentioned above, and ( cmm), the average 02-02 interaction (i.e., the lateral interaction). Plots of these quantities (in reduced units obtained by dividing energy by cos, the site-site interaction with the solid, are shown in Figures 2 and 3 as a function of temperature. Values of the total average potential energy per molecule ( U ) = (eps) + (e-) are listed in Table I1 together with relevant data for each of the simulation runs. The most striking feature of the plots is the difference between the patch and the full surface simulations. In this study, lI6of the atoms are on the edges of the patches at 0 K and this fraction becomes larger as the periphery of the patch becomes nonlinear at finite temperature. Even more important at finite temperature is the fact that the patches can expand whereas the full surface is constrained to a fixed density except at very high temperatures where molecules can leave the monolayer. Thus, the rapid decrease of the (negative) lateral interaction energy with increasing temperature for the patch systems is due to this expansion. On the other hand, the differences in average
584 Langmuir, Vol. 3, No. 4 , 1987
Bhethanabotla and Steele
10.0* 20.3* 25.0 29.0 32.0 34.0* 35.0 38.0
Table 11. Simulated Systemsa time steps after equilibra( U )/e,,tion temp ( U )/eoR Full Surface -43.60 2500 40.0 -40.66 2000 42.1* -40.38 -42.92 44.6* -40.17 9000 -42.53 -42.29 16000 45.0 -40.08 11000 55.3* -39.13 -41.97 60.0 -38.69 3500 -41.85 -37.65 10000 75.0 -41.69 -40.85 10000
20.0* 25.6* 28.6* 30.0 31.7-
-43.38 -41.86 -41.13 -41.22 -40.14
temp
.."",
time steps after equilibration
35.0 36.2* 39.4* 60.0 70.0
-40.25 -39.82 -39.31 -36.67 -35.65
L
I
Patch
13250 3000 3000 10250 3000 5000 5000
Patch 2000 4000 3000 20000 3000
I
10000 3000 3000 20000 10000
0
0.00,
I.,
1
.I
I
r
Full surface
A
n
30
45 Temperature
60
Figure 4. Translational order parameter, which is the structure factor evaluated by using the 0 K inverse lattice parameters, is shown here as a function of temperature. Thermal expansion causes a gradual decrease in OP2, but the sharp drops are associated with melting.
"uns at constant N,V,E are denoted by asterisks; others are at constant N.V,T. Patch
molecule-solid energy observed at high temperature were unexpected. However, the energy minimization calculat i o n ~and ~ , the ~ simulated orientational order parameters discussed below show clearly that there is a strong preference for surface-parallel orientations in these systems, and molecules in an expanding patch have less difficulty in maintaining this parallel orientation at high temperature than the crowded molecules in a full-surface monolayer. To locate phase transitions in these systems, it is helpful to evaluate the appropriate order parameters. The parameters calculated here have been defined and evaluated in the previous N,/gcb simulation^.^^^ One has two orientational order parameters denoted OP, and OP, for the out-of-plane and in-plane ordering respectively. These are OP4 = (cos 2 a )
---O --I-0-
0
-
gn 0.50 I
o.oo,5
30
45 Temperature
\]
I 1
(Y2 cos2 p - 7 2 )
1 I
I
where /3 is the angle between an axial vector and the normal to the surface. In this case, all molecules become coplanar with the surface at 0 K so that the limiting value of OP, is -1/2. Translational order is often characterized by evaluating the structure factor S(k). For a two-dimensional system,
S(k) = (eck.rfi) where T , ] is the separation vector between the centers of mass of molecules i and j . In this study, we evaluated S(k) only for the k vectors that characterized the periodicity of the original (low-temperature) phase. Thus, we took k* = (0.75, -0.3079) and 0, 0.6158), with reduced values defined as h,* = k, (2.46 A), i = x or y. The structure factor for these values of k is denoted by OPz. The temperature dependence of these simulated order parameters is shown in Figures 4-6 for both the full surface and the patch. First, note that translational order is lost abruptly at T = 27 ( f l )and 36 ( f l )K for the patch and the full surface. Of course, as the patch expands with increasing temperature, OP2 will decay even if the system remains solid; however, this decay will be gradual and indeed, the difference between the low-temperature slopes of OP2 for full surface and patch is a consequence of the
60
Figure 5. Out-of-plane orientational order parameter is shown here. The gradual decreases with increasing temperature are due to an increase in the amplitude of librational motion, but the sharp change and subsequent patch-full surface difference are connected with melting in the full-surface system.
where a is the angle between the in-plane projection of a molecular axial vector and the X-axis of the simulation box. Since all molecules align with the Y-axis in the limit of 0 K, the limiting value of OP, is -1; also,
OP3 =
Full surface
g*o.501 I
O.0Ol5
0
Patch F ~ I Isurface
45 Temperature
30
60
Figure 6. These plots of the temperature dependence of the in-plane orientational order parameter show that both solids remain strongly ordered u p to their melting points. However, significant ordering remains in the liquids just above their melting points. OP, is defined to be (cos 26a) where ba is the angle between projections of the molecular axes and an in-plane vector lying parallel to the sides of the graphite hexagons (and thus, perpendicular to the molecular axes in the initially aligned O2 crystal). Nonzero values of this parameter in the liquid thus indicate a preference for the molecular axes to lie perpendicular to the sides of a graphite hexagon in a direction correlated with their initial orientations.
expansion. The evidence of Figure 4, taken together with the trajectory plots shown below, provides convincing proof that melting to a two-dimensional liquid is occurring. It is gratifying to note that the patch melting temperature
Langmuir, Vol. 3, No. 4, 1987 585
Simulations of Oxygen Monolayers on Graphite
I
Patch,
I
30K
Ah A A
A
a
d
A
. A
1
Full surface, 38K
d
A A
A
I
1
I
A
32K 0
A
.
I
A
a
-
I
A
2.51. A
A
A
h
I
cos p
Figure 7. Distribution functions for the out-of-planeorientation of the molecules in the full surface system are shown here. Plots are for cos p, where /3 is the angle between the surface normal and the molecular axes. Although the peaks show a strong preference for orientations parallel to the surface, it is evident that a significant fraction of molecules are standing on end at 38 K. (The.experimenta1evidence is that further increases in the density will give rise to new ordered phases which consist entirely of molecules standing on end.5)
agrees nicely with the experimental value of 26 K (see Figure 1). Experiments also show an increase in the melting temperature for the full layer, which first occurs a t a coverage =7.4 X molecules/A2. However, the estimated experimental melting point is not quite as high as that indicated by the simulation. In sharp contrast to the behavior of Nzmonolayer on graphite, the plots of the orientational order parameters in Figures 5 and 6 show that the simulated O2 layers remain orientationally ordered up to the melting points. However, some in-plane ordering appears to be present in the liquids a t temperatures just above the melting point and, surprisingly, the out-of-plane order parameter for the filled surface undergoes a jump a t or near to the melting point. Evidently, this jump corresponds to that shown in the average molecule-solid energy plotted in Figure 2. This appears to be the first example of a two-dimensional disordering transition that simultaneously involves a sharp change in an order parameter for the third dimension. It is interesting to note that Fain and co-workers5have observed that the melted &phase possesses considerable inplane order, especially for the full surface systems. They suggest a distinct &phase, as shown in Figure 1and have also speculated that changes in the out-of-plane order may be required to explain their observations. Of course, a much more detailed picture of the orientational ordering can be extracted from the simulations by calculating the full distribution function for the in-plane angle a and the out-of-plane variable cos p. Distributions of cos @ are shown in Figure 7 for the full surface a t a temperature below (32 K) and above (38 K) the melting transition. Aside from a general thermal broadening, it is evident that a significant fraction of molecules in the high-temperature phase are nearly on end with cos (3 = f l and that this behavior is absent at the lower temperature. A few distributions of in-plane angle are shown in Figure 8. In particular, the top two curves for the liquids just above their melting points exhibit significant orientational ordering in the preferred direction for the solid phase.
"k
586 Langmuir, Vol. 3, No. 4, 1987 w
Bhethanabotla and Steele
Y
a2
I !. r.i.
I
.......... ............. ............I ............. i ........... 1 ............ *I ........... 1............ ............I I. ...........I ...... >\...., ............ i.....**& ..., ~
c
L .
.p.i4?i
tirnl
;.,---*...,.{ ,.LI.L.J',\J
(P,CO..d
Figure 11. Mean-square in-planedisplaeementsare plotted here for the full-surfacesystem. Reduced units are defied by dividing distance by a = 2.46 A The linear variations with time that appear to characterize the high-temperature curves indicate diffusional motion in the liquids. The slow rises in the lowtemperature curves for the solids are believed to be due to the floating of the crystals across the surface.
1
I
Figure 9. Computer-generated plots of the in-plane center-ofrima trajectories me shown here for the 192 molecules in the patch system at three temperatures, one below the melting point and two above. The total time span of these trajectories is -50 ps. The results clearly show twedimensional liquid-vapor equilibrium at 60 K as well as a highly ordered solid at 25 K. The use of periodic boundary conditions means that these pictures are part of an infinite array generatedby translating the central box along the x and y directions. 1 .
"
.,.
,
-2 On
I 2 lima (pisosr)
3
Figure 12. Logarithmic plots of the reorientational correlation for the full-surfacesystem. The function plotted is (cos am@)), where W t )is the change in the in-plane axial orientation angle in time t. Linear decays on this scale are characteristic of rotational diffusion, hut the nonlinear regions indicate that memory effecta are significant in these systems. Restricted libration gives rise to the nondecaying curves for the solids at 25 and 35 K.
Figure 10. Trajectory plots for the in-plane motion of the molecules in the full-surface system are shown here. (As in Figure 9, periodic boundary conditions are used.) The floating motion of the solid at 32.0 K is clearly shown. Here the displacements occur along the diagonal of the rectangular unit cell. The liquid at 38.0 K shows a Considerable amount of translational order in the central region of the picture. Conceivably,this could be a domain in a two-phase equilibrium state, but the small size of the computer sample precludes any definitive argument along this line. Even at 45 K, remnants of preferential displacements to neighboring unit cells can be seen. the surface at the edge of an absorbent crystallite. Just above the melting point at a temperature of 38 K, the trajectories indicate that considerable order remains in the fluid phase, and some remnants remain even at 45 K, thus supporting the suggestion of an unusually ordered fluid in this part of the phase diagram. Although the time span for the trajectories is reasonably long (50 ps) on a molecular scale, it is likely that the order would be washed out if observations were taken over a much longer intervalhowever, this by no means precludes its measurement in a suitable experiment. The velocity time-correlation functions for translation perpendicular to the surface in these monolayers indicate that this degree of freedom is nearly harmonic oscillation.
However, the nature of the motion parallel to the surface is much more variable, changing from oscillatory in the solids at low temperature to diffusional in the high-temperature fluid phases. To illustrate this point, mean-square displacements are plotted in Figure 11for the fully covered surface. Only the motion parallel to the surface is shown, since perpendicular displacements are invariably small and independent of time after an initial rise, as expected for harmonic oscillation. The onset of diffusion parallel to the surface at a temperature between 35 and 45 K is clearly indicated in Figure 11. The slowly increasing displacements shown for the solids at 25 and 35 K are due to the floating motion of the entire layer across the surface. One of several ways of characterizing orientational dynamics in these systems is to evaluate the time-dependent changes in the direction cosines for the molecular axes. Here again the in-plane and out-of-plane dynamics are very different, with the out-of-plane reorientation consisting of a nearly harmonic libration for all systems simulated. In-plane, one evaluates (cos - c i ( O ) ) ) , where CY is, as above, the orientation angle for in-plane projection of the axial vector. In an orientationally ordered solid, this time-correlation function will quickly reach a constant value close to unity, corresponding to small librational displacements. On the other hand, the function will decay exponentially to zero if rotational diffusion is present:' as one might expect for the orientationally disordered fluid (21) Steele, W.A. Ado. Chem. Phys. 1976,34, Sed. IIIB.
,
Simulations of Oxygen Monolayers on Graphite phases. The logarithmic plots of this time-correlation function shown in Figure 12 provide added confirmation to the general picture developed here. Below the melting point for the full surface, the monolayer is orientationally ordered and the molecules are undergoing librational oscillations in the plane of the surface. Above the melting point, the motion is reasonably well approximated by rotational diffusion with a rotational diffusion constant (proportional to the slopes of these plots) which increases with increasing temperature.
Discussion In general, the simulations reported here are in agreement with the experimental findings for the 02/gcb system. Melting points for the partly covered and the fully covered surface are close to those observed, although some discrepancy between simulation and experiment is observed for the full-surface melting point. Note that a very strong temperature variation of melting has also been observed for the nearly commensurate layers of Nz on gcb which is reproduced by ~imulation.~ It is thus somewhat surprising that the experimental dependence of melting point upon density is so small for Oz. The melting process in the Ozpatch system is reasonably straightforward-it is first order, with a finite energy and density change associated with the transition. From the plots in Figure 3, the energy change can be roughly estimated to be 0.6 reduced units or 190 J/mol (experimental values range from 170 to 275 J / m 0 1 ) . ~It~ is not feasible to estimate the density change because of problems in locating the liquid-vapor dividing line. Melting in the filled surface system is more interesting. As is well-known from standard bulk thermodynamics, the first-order isothermal phase changes that occur in isobaric systems are converted to second-order nonisothermal processes in isochoric systems. At constant volume, one observes only a change in slope of the energy vs. temperature plot (i.e., a discontinuity in heat capacity) which is due to the gradual liberation of the energy of melting with increasing temperature. In our quasi-two-dimensional case, a constant areal density is maintained so one might reasonably expect a two-dimensional analogue to the three-dimensional isochoric case. Certainly, no discontinuity is observed in the lateral interaction energy E,, a t the melting point. However, a change in the third dimension of the configurations is available to physisorbed systems-for example, molecules could be promoted from the first to the second layer or, as in the present case, molecular orientation relative to the surface can change. The shift to more upright orientations gives more room for the molecules to become liquid and appears to produce an isothermal (or nearly so) melting process with a finite energy change, estimated from the curves in Figure 2 to be 0.55 reduced units or 175 J/mol, with a rather large uncertainty. To our knowledge, this is the first time that such a process has been observed. It is gratifying that the simulations show successive solid, liquid, and vapor phases as temperature increases. There is some controversy in the literature concerning the value of the two-dimensional liquid-vapor critical temperature for O2 on gcb, with 60 and 65 K having been ~ u g g e s t e d . ~ A precise determination by simulation would be extremely difficult and thus has not been attempted. The trajectory
Langmuir, Vol. 3, No. 4, 1987 581 plots indicate a value between 60 and 70 K, with some weak indication that an intermediate temperature rather than one close to 60 K may be correct. However, it must be remembered that all values of the temperature quoted in this paper are dependent upon the choice of well-depth parameter for the 02-O2 interaction. There are experirnentall4 and theoreticalzzindications that this well depth may be significantly reduced from the bulk-phase value used here due to substrate-mediated effects. If so, the temperature scale for the phenomena dependent upon lateral interactions should be scaled. To a first approximation, the scaling is the same (but inverted) for the well depth, for the temperatures of melting, critical point, etc., and for the energies of lateral interaction. The molecule-solid interaction remains unaffected by any such changes, which is pleasing in light of the excellent agreement between simulated and experimental properties that are dependent upon this potential function. One of the most striking difference between O2 and Nz adsorbed on gcb is in the melting points of the patch (or partly filled) surface. The solid is strongly stabilized for N2 which has a melting point almost double that of O2 but a site-site well depth that is N * / less. ~ This behavior is due to a combination of the stabilization energy due to registry over the graphite lattice and to the presence of a strong quadrupole-quadrupole interaction in the N2. In fact, it can be argued that the N2 melting point is higher than the two-dimensional liquid-vapor critical temperat ~ r eso, that ~ ~ liquid is not observed a t any temperature for these layer densities. In the O2 case, it is reasonable to conclude that the formation of liquid at extremely low temperature (relative to c / k for the 02-02 interaction) gives rise to strong and long-range positional and bondorientational order. Although we found little evidence for a discrete &phase such as that suggested in Figure 1, it is clear that the order in these liquids at 30 K is unusually strong. Floating motion of the two-dimensional crystal has been observed previously4 in a simulation of the unaxially compressed N2 layer. It is interesting to note that floating occurs only in a specific direction for the O2 and N2layers. Both layers are incommensurate but still possess angular order relative to the underlying graphite lattice. Intuition and theoryz0suggest that floating must occur along a direction of constant (or almost constant) potential energy. Given the detailed interaction law used in these simulations, it would be of interest to calculate the energies of the N2 and O2 crystallites as a function of position relative to the underlying lattice.
-
Acknowledgment. Support for this work was provided by Grant DMR84-19261 from the Division of Materials Research of NSF. Helpful conversations with S. Fain, R. Etters, M. Cole, and M. H. W. Chan are also acknowledged. Registry No. 02,7782-44-7; graphite, 7782-42-5. (22) Sinanoglu,0.; Pitzer, K. S. J. Chern. Phys. 1964,41,1322. Bruch, L. J . Chern. Phys. 1983, 79,3148. McLachlan, A. D.Mol. Phys. 1964, 7, 381. (23) Migone, A. D.; Chan, M. H. W.; Nilscanen, K. J.; Griffiths, R. B. J. Phys. (Park) 1983, C16,L1115. Butler, D. M.; Litzinger, J. A.; Stewart, G. A.; Griffiths, R. B. Phys. Reu. Lett. 1979, 42, 1289. (24) Marx, R.; Christoffer, B., preprint. See also ref6, Stoltenberg and Vilches.
