Simulation studies of bilayers of N2 adsorbed on graphite at 25 and 35 K

Surface Science 171 (1986) 83-102 North-Holland, Amsterdam

SIMULATION STUDIES OF BILAYERS GRAPHITE AT 25 AND 35 K Alexei VERNOV

83

OF N 2 ADSORBED

ON

* and William STEELE

Department of Chemistry, The Pennsylvania State University, Universi(v Park, Pennsylvania 16801, USA Received 30 August 1985; accepted for publication 2 December 1985

Molecular dynamics simulations of bilayers of model nitrogen molecules adsorbed on the graphite basal plane at 25 and 35 K are reported. Systems considered include commensurate bilayers and three uniaxially compressed films, namely, two with different compressions along the X-axis (perpendicular to the herringbone lattice glideline) and one with a compression along the Y-axis, relative to the commensurate spacing. For comparison, commensurate and X-compressed monolayer films were also simulated. Properties simulated include: center-of-mass (COM) molecule surface densities as a function of separation distance; in-plane COM pair correlation functions; in-plane and out-of-plane distributions of the molecular axis orientations; average potential energies for Nz-surface and for N2-N 2 interactions; and velocity time-correlation functions for translational and orientational motion perpendicular to and parallel to the surface. These results indicate that the N 2 layers at these temperatures are nearly harmonic oscillator solids, with considerable orientational freedom. Molecules tend to be coplanar with the surface in both layers. The in-plane orientations for the compressed films form a herringbone lattice, in both layers, with variable angles between the molecular axis and the crystal glidelines; the commensurate film does not show herringbone ordering for the potential model used in this work. It is concluded that both the orientational and the translational degrees of freedom are strongly coupled, between layers as well as within the layer. It is also shown that the Y-compressed film is unstable and that both layers in the bilayer spontaneously rearrange into a (rotated) uniaxially X-compressed structure.

1. Introduction Studies of the adsorption of nitrogen on graphitized carbon blacks are quite e x t e n s i v e , i n p a r t b e c a u s e n i t r o g e n is s o w i d e - s p r e a d i n its u s e as a n a d s o r b a t e in c o n n e c t i o n w i t h B E T s u r f a c e a r e a m e a s u r e m e n t s , b u t a l s o b e c a u s e n i t r o g e n is a p r o t o t y p e s i m p l e n o n - s p h e r i c a l m o l e c u l e w h i c h h a p p e n s t o e x h i b i t i m p r e s s i v e l y c o m p l e x p h a s e b e h a v i o r f o r its m o n o l a y e r a d s o r p t i o n o n g r a p h i t e [1-21]. At low temperature, several two-dimensional solid structures are now * Permanent address: Department of Physical and Colloid Chemistry, People's Friendship University, Ordzhonikidze St. 3, 117302 Moscow, USSR.

0039-6028/86/$03.50 © E l s e v i e r S c i e n c e P u b l i s h e r s B.V. (North-Holland Physics Publishing Division)

84

A. Vernov, W. Steele / Simulation of N 2 adsorbed on graphite

known, including a commensurate ( f 3 × ~ - R 3 0 °) and two incommensurate solids, all of which transform from orientationally ordered at low temperatures to disordered phases at temperatures near 30 K. A number of experimental techniques have been employed to study this system, especially at low temperatures where conventional isotherm measurements are not feasible. Among the most productive of these are LEED [3,6-8] and heat capacity measurements [2,10 12]. As a supplement to the experiments, computer simulations give one the possibility of determining microscopic ordering and molecular distributions, if a realistic model for the intermolecular interactions can be devised. Indeed, several such models have been employed in rather extensive simulations of the N 2 graphite system [14 17,20]. These models differ in their detailed parameterization but have in common a molecule solid energy which reflects both the non-spherical shape of N 2 and the periodic nature of the graphitic surface. In addition, the N 2 - N 2 interaction function is chosen to be one which is known to give a good account of the properties of bulk solid and liquid N 2. Simulation studies already reported include the commensurate monolayer at temperatures from 5 to 45 K [14,17]; unaxially compressed layers at 15 and 40 K [20]; and the fluid at sub-monolayer and multilayer coverages with T = 74 K [16,21]. One of the findings of the last of these was that second- and third-layer molecules have a significant effect upon the density and orientation of the first-layer molecules (which obviously leads to energetic and dynamical property changes in these molecules as well). Inasmuch as the compressed monolayer phases at lower temperatures show complex changes in orientational ordering, it seemed that a comparative study of monolayer and bilayer films at low temperature might be of interest. Such simulations are reported here. Temperatures of 25 and 35 K were chosen, and uniaxially compressed bilayers (parallel to and perpendicular to the glidelines of the herringboneordered phase) were simulated, together with the commensurate bilayer and, for completeness and for comparison with previous work, some relevant monolayer systems. It will be shown that the presence of the second layer does indeed affect the first-layer properties of these systems and furthermore, that molecular orientations and positions in the plane parallel to the surface are strongly coupled in these films.

2. Description of model and simulations As in previous simulations of the N2/graphite system, the intermolecular interactions were taken to be the pair-wise sums of site site Lennard-Jones 12-6 functions plus a electrostatic quadrupolar energy based on three discrete charges located at the N atoms and at the center of symmetry (with opposite sign and twice the magnitude of the atomic charges). Values of the well depths

A. Vernov, W. Steele / Simulation of N 2 adsorbed on graphite

85

Table 1 Parameters of the model potential energy ONC (nm)

fNc/k (deg) l (nm)

0.3359 31.9 0.1110

ONN (nm)

I[NN//k (deg) q*

0.3318 36.4 18.676

q* = q,/cV/~uca, a = 0.246 nm.

(E/k), site diameters (o) and the atomic charge used here are listed in table 1. This model for the molecule-surface energy involves a pair-wise sum over the carbon atoms of the graphite substrate. It is convenient to express this as a Fourier series which reflects the periodicity of the surface lattice [22]. As is known from previous studies of this system [23], only a few terms need be included to obtain an adequate representation of the periodic part of the molecule-solid potential and also, relatively simple expressions can be written for the leading (non-periodic) term in the series. Although there are significant indications that the sum over a t o m - a t o m functions does not correctly model the anisotropy in the graphite planes of the substrate and thus yields a periodic component of the potential which is too small (possibly by as much as a factor of two) [24], our choice of model was determined by the fact that comparison with previous simulations of this system is only possible if the a t o m - a t o m sums are employed. Another difficulty with the model arises when choosing the magnitude of the partial charge q or equivalently, the N 2 quadrupole moment (ql2/2, where l is the site-site separalion). Subsequent to the first simulation study in this series, this charge was reduced by - 30% to give an "effective" quadrupole moment which disagrees with the experimental measurement at low density but which appears to give the best property values for the dense bulk phases of nitrogen [25]. Evidently, it is impossible to choose q to be identical with that used in all previous studies, so the most recent value was taken. (It will be shown that the change has a marked effect on the orientational disordering temperature for the commensurate monolayer.) Finally, no attempt was made to include three-body energies, particularly those involving substrate-mediated interactions. The justification for this is simply that their inclusion would make the simulation to be impossibly slow. Since the N 2 - N 2 potential used is an empirical representation of the effective interaction in the dense phase, it is not clear how one should include the substrate-mediated modification in any case. Table 2 gives a listing of the simulation runs completed and analyzed in this work. Because of the large energy jumps associated with interlayer transfer of an adsorbed molecule, a constant temperature molecular dynamics algorithm was used, based on the suggestion of Evans et al. [26] for atomic systems; minor modifications were necessary to deal with the rotational kinetic energies of the N z molecules. However, the shift to constant kinetic energy from the

16 16

~} A f t e r equilibration. b} C o m p r . = Compression. c) C = C o m m e n s u r a t e .

18 18

8 9

16 16 16 16 16 16

16

18

18 18 17 17 16 16

1

/~, ( 2 / f 3 )

1~

2 3 4 5 6 7

No.

1.063 1.051

1.000 1.000 1.059 1.059 1.125 1.125

1.000

Nf

Table 2 Description a n d results of simulation runs

0.938 0.949

1.000 0 1.059 1.059 0 1.125

1.000

N~

Y-uni Y-uni

C ° C ° X-uni X-uni X-uni X-uni

C c}

C o m p r . b)

25 35

35 25 25 35 25 25

25

T (K)

11400 15425

13300 10575 14450 10950 16775 10600

13425

Timesteps ~I

- 35.97 35.34

-36.01 -36.19 - 36.21 - 35.57 -34.73 - 35.78

- 36.48

U~*~,1

- 5.15 -4.90

- 4.86

- 5.16 4.92

-5.16

- 5.37

b*~,2

- 9.04 8.76

-8.82 -9.07 - 9.20 - 8.93 -9.39 - 8.65

- 9.09

Ul*l

-8.88 -8.53

-10.35

-9.85 - 9.70

-9.33 -9.13

U2"2

Q

g

A. Vernov, W. Steele / Simulation of N, adsorbed on graphite

87

usual constant energy M D turned out to be unnecessary, since only one molecule underwent interlayer transfer in the entire study (in run No. 9). In fact, this is a significant result: namely, at the low temperatures considered here, the probability that a molecule will change from one layer to another is negligibly small on a picosecond timescale. (This is in contrast to the situation at 75 K.) Each system was equilibrated for 3000 timesteps, with a timestep in the predictor-corrector algorithm of length 2.02 × 10 -3 ps. The data were then gathered over runs of varying duration, as noted in table 2. Averages obtained over blocks of 500 timesteps indicated that orientational equilibrium was reached rather quickly. Although the attainment of translational equilibrium within a layer is more difficult, the results indicated that this was also achieved - see particularly the discussion of the Y-compressed uniaxial phase given below. Reduced units for the energies given in table 2 are energy/cNc (per particle); reduced units for the distances reported below are distance/0.246 nm, Initial configurations for all systems were the herringbone lattice of N 2 molecules coplanar with the surface and with if, the angle between the molecular axes and the crystal glideline, equal to 45 °. However, the secondlayer lattices were translated relative to the first layer such that each secondlayer molecule was initially located over the center of a triangle of first-layer molecules. First- and second-layer molecules were initially set at z = 0.332 nm and 0. 676 nm respectively. As in previous simulations, periodic boundary conditions were used in the x- and y-directions and molecular velocities in the direction were instructed to reverse at a plane distant from the surface however, molecules did not actually leave the surface phase in these systems. Both time correlation functions (TCFs) and equilibrium distributions of molecular C O M position and orientation were calculated. The time correlation functions were all single-molecule quantities and were averaged over multiple time origins as well as over all molecules. The equilibrium singlet and pair distributions to be reported below were also averaged over time as well as the molecular configurations at a given time.

3. Results 3.1. Translational and orientational order

At the low temperatures of these simulations, the C O M equilibrium positions and dynamical properties are those of a nearly harmonic librator-oscillator solid, as was shown for monolayers in the previous study [14,27]. The distance dependence of the molecule-solid singlet density shows the expected separation into layers, with a small mean square amplitude of oscillation

88

A. Vernou, 14A Steele / Simulation of N 2 adsorbed on graphite

0.3

#6

#1

#8

0.2 v

0.1

1.~ I

2 z ~

3

0

I

I

2

3

I

0

I

2

3

4

z~

Fig. 1. Center of mass (COM) molecular densities as a function of distance (in reduced units) from the graphite surface. Numbers indicate the simulation run according to the code in table 2.

perpendicular to the surface. Selected curves for this density at T = 25 K are shown in fig. 1. In fig. 2, the velocity TCFs for motion perpendicular to the surface are shown. They appear to correspond to slightly damped oscillators. However, both the frequency and the damping are significantly affected by the presence of a second layer and, not surprisingly, these constants for secondlayer motion are quite distinct from those for the first. It is notable that the frequency of oscillation for a first-layer molecule against the surface goes up (as expected) when a second layer is imposed, but the damping also increases. T C F s for the torsional motion associated with rotation of the molecule axes in the plane perpendicular to the surface are shown in fig. 3. Here again this motion for first-layer molecules is dominated by the large changes in molecule-solid potential energy, but the modification of the first-layer T C F due to the presence of a second layer is quite significant. The strong damping of the second-layer T C F is characteristic of motion in a bulk solid and may be interpreted as the result of a wide distribution of librational frequencies. Inasmuch as the translational velocity T C F Cv(t ) has the appearance of harmonic oscillation, an attempt was made to fit the well-known result for a Brownian harmonic oscillator curve to the simulation. This failed, in part because of anharmonicity in the potential (which is a sum of inverse 10th and 4th powers of the distance and thus, not very harmonic for large displacements), but mostly because the damping was not correctly given by the simple Brownian motion model. It appears likely that memory effects in the random part of the forces are significant on the timescale of these calculations. Some plots of the in-plane velocity TCFs are shown in figs. 4 and 5 for the same systems as in figs. 2 and 3; namely, commensurate layers at 25 K. Here, the changes associated with changing coverage are much less noticeable and

A. Vernov, W. Steele / Simulation of N. adsorbed on graphite

1.0

I

I

Cv(t'0.5t ) (perp.) 0 -0.5 -I.0!

89

I

t

4.\"

I

-

IJ~..(i/1) I 0.5

I 1.0

time (picosec)

1.I 5

2.0

Fig, 2. Time correlation functions for the component of the translational velocity perpendicular to the surface. All curves are for commensurate layers at 25 K. The dashed curve shows the TCF for a monolayer (1/1) and the two solid curves are for first-layer molecules (1/2) and second-layer molecules (2/2) in the bilayer (run 3 and 1 in table 2, respectively).

the decays to zero are more pronounced, as expected for motions which are primarily controlled by the N 2 - N 2 part of the potential. Another time-dependent quantity evaluated in these simulations was the mean square molecular displacement, in directions in the plane and perpendicular to the plane of the surface. There was clear indication of diffusional motion parallel to the surface for the simulated commensurate and bilayer films at 74 K, thus throwing doubt on the suggestion that solid layers could be present at this temperature. We omit the plots, but all mean square displacements evaluated here reached an essentially constant value after 1 ps, at 25 and 35 K, indicating that these layers are indeed (nearly) harmonic oscillator solids. Ample confirmation of this finding is obtained from the equilibrium pair-wise C O M distribution functions. Several such curves are given in fig. 7. The lattice spacings and coordination numbers for the perfect commensurate lattice (part a) and for the X-compressed lattice (part b) are shown by vertical lines, with heights proportional to the number of neighbors to a given molecule at distance r divided by r. It is evident that both the commensurate

90

A. Vernov, W. Steele / Simulation of N_, adsorbed on graphite

I'0~!t I1%

I

I

I

1.0 time (picosec)

1.5

I~\ f(z/2) I li~'/('") I Ilk"/(~lZ)

X

f,

0

0.5

2.0

Fig. 3. Angular velocity TCFs for libration in the plane perpendicular to the surface. Curves are shown for the same systems as those of fig. 2.

I.O~

I

I

I

h

0.5~

Cv(t)

olk

(parallel)I

-0 5 ~ • 0

I

0.5

(112)

(212)

1.0 time (picosec)

1.5

2.0

Fig. 4. Translational velocity TCFs for motion parallel to the surface in the same systems as those of fig. 2.

A. Vernov, W. Steele / Simulation of N? adsorbed on graphite

1.0

I

I

91

I

0.5 C(u(t) (parallel)

L

0

ll/t)

(2•2) -0.5 -

0

I

0.5

I

1.0 time (picosec)

I

1.5

2.0

Fig. 5. Angular velocity TCFs for reorientation in the plane of the surface for the systems of fig. 2.

and the X-compressed monolayers have structures that are very close to those for the perfect lattices, and exhibit only a moderate degree of thermal broadening. In part c of fig, 7, the Y-compressed monolayer is compared with an X-compressed layer with essentially the same density. Although the Y-compressed layer is less perfect, as shown by the increased peak broadening, there are strong indications that the two compressions lead to lattice that are very similar. In fact, the weight of the evidence in this study leads to the conclusion that this Y-compression yields an unstable system which rearranges into a rotated X-compressed lattice, at 25 and 35K. This point is illustrated convincingly in the two snapshots of the Y-compressed first layer shown in fig. 8. One is of the initial configuration, with glidelines parallel to the Y-axis, as shown by the dashed line. The second snapshot was taken at the end of the simulation and shows quite clearly that the glidelines have rotated. Furthermore, the distances between molecules along the x- and y-axes are tabulated in table 3 for the initial and final states of the simulation and show that not only have the molecules reoriented into a new crystal, but the spacings have changed so that the final state is indeed a single-crystal X-compressed layer. In a previous simulation of uniaxially compressed films [20], a similar rearrangement of the Y-compressed layer was observed, but the final state appeared to be polycrystalline X-compressed with significant domain wall structure. The chief difference is that the simulation runs were considerably

92

A. Vernou, I~ Steele / Simulation of N 2 adsorbed on graphite

Cv(

0

0.5

I.O

1.5

time (picosec) Fig. 6. Translational velocity TCFs for the X-compressed monolayer (:~6 in table 2). All three cartesian velocity components are given in order to show the anisotropy in velocity produced by the compression.

longer in the present study, which indicates that the previous result may have been for an intermediate state during the transition from Y-compressed to X-compressed. An interesting contrast to these results is given by the translational velocity TCFs shown in fig. 6 for the uniaxially compressed monolayer at 25 K. Motion in both the z- and y-directions is hardly changed from the corresponding commensurate layer, but the X-compression has produced considerably more oscillatory character for motion in that direction. The angular velocity TCFs shown in fig. 9 for this system are less influenced by the uniaxial compression, but the in-plane rotation around the z-axis is considerably more oscillatory than the corresponding rotation in the commensurate monolayer (see the dashed curve in fig. 5). The question of the nature of the in-plane orientation ordering in the first and second layers of the simulated systems can be most easily answered by inspection of pi(d~), the probability density that a molecule in the ith layer has angle ~ measured relative to the x - y coordinate axes. A number of these curves are shown in fig. 10. The set in the left-hand column are all for 25 K

A. Vernov, W. Steele /Simulation of N, adsorbed on graphite I

I

!

I

I

1

I

I

I

#1

i

I

I

I

I

-#6

a

8O

2

I

I

I

-¢#4 ---4#'8

b

r ~

93

c

4 r~

80

2

-

4 r~

8

Fig. 7. Center of mass pair correlation functions projected into the plane of the surface for first-layer molecules in the simulation runs at 25 K denoted by the numbers. The vertical lines show the separation distances in the perfect two-dimensional crystals; the heights are proportional to (coordination n u m b e r s ) / r * . (a) Run 1 is a commensurate monolayer, (b) run 6 is an X-compressed monolayer, (c) run 4 is for the first layer in an X-compressed bilayer and run 8 (dashed curve) is for the first layer in an Y-compressed bilayer.

and reveal that the commensurate layer is already above the herringbone disordering transition, since one sees three weak peaks (in both the first- and second-layer curves) due to a residual ordering in the cells of six-fold symmetry. In fact, there is a mild tendency for the molecular axes to point toward their nearest neighbors in this system. A slight compression in the x-direction (system ~ 4 ) causes pronounced herringbone ordering to appear in both layers of the bilayer system and in the monolayer system (dashed curve). The most probable angles and the degree of ordering change as the X-compression increases, as shown by system :#:6. The ordering in the second layer is less pronounced than in the first but it is interesting to see that the herringbone ordering persists in this system to T > 35 K. The curves of p(~) shown for the Table 3 Spacings in the Y-compressed monolayer

6X*/3

Y*/~/-3-

Initial

Final - relative to new glideline

1.000

0.960

0.942

0.990

94

A. Vernou, W. Steele / Simulation

/ / /- ,

/

4

OO

-

Initial --

8

/

\ 6--\

/ ,

\

’ \ ’ \

/

1 1

‘I

2

t

I

Oo

‘,‘/I \

/

\

,I’

\

,

/ / I\,\,\,‘-

‘!, /’

4’

/I

/ ’

I

/

I\

/I

/ I\

; \I p 8

original glideline

#I

State

/

‘I

l\

I ’

f

16

I’,‘,’

4_\,\ \

I

12

6

of N, adsorbed on gruphrte

’ \

\

/ /

/

’ \

,

’ ‘\

/

\

/

\

/‘/‘/’ 1

\

\

/

, I2

\

\

/

\

I\

i

1

I6

X*

Fig. 8. Two vertical views of the instantaneous positions and orientations of N, molecular axes in the computer box are shown. Only first-layer molecules are pictured; the initial snapshot shows the Y-compressed monolayer and the final snapshot (after = 13000 timesteps) shows the spontaneous reorientation and translation to an X-compressed lattice with a rotated set of glidelines.

Y-compressed system at the bottom of the array can only be understood if it is realized that the double peaks spaced by roughly 90” correspond to a herringbone structure, but ordered in a frame with glideline originally at + = 90”, but finally rotated by 60” to a 30” angle relative to the original axes. For the first layer, these results provide solid confirmation of the argument that the Y-compressed system is unstable relative to an X-compressed herringbone-ordered lattice. The ordering is not as sharp in the second layer, but the double peaks in this case are centered about the same angle and split by the same amount as in the first layer. A comparison of the sets of curves for 3.5 K shows that the Y-compressed system shows only weak herringbone ordering

A. Vernov, W. Steele /Simulation of N 2 adsorbed on graphite

1.0

I

I

I

1.0

1.5

95

I

0.5

I

l

U II II -0.5

0

0.5

time (picosec) Fig. 9. Same as fig. 6, but for the angular velocity TCFs. corresponding roughly to the (rotated) X-compressed structure in this case. It is likely that domains and disordered regions are present in this case. Fig. 11 shows the out-of-plane ordering at 25 K for the commensurate and the two X-compressed bilayer systems. For the first layer, a comparison between bilayer and monolayer distributions of cos fl is given. Here, fl is the angle between the molecular axis and the surface plane. It is evident that these adsorbed molecules tend to lie parallel to the surface in the second layer as well as in the first. Also, the presence of a second layer has the effect of flattening the molecules below and effectively eliminates those molecules which might be standing on end. The ordering in the second layers decreases significantly from that for the commensurate system as more molecules are added, but the flattening effect on the sublayer remains. It is reasonable to conclude that this increase in coplanarity restricts the in-plane orientations and thus causes the sharpening of the q~-distributions shown in fig. 10. 3.2. A v e r a g e energies a n d in-plane angles

It is now possible to understand the differences in the average energies reported in table 2. The first question is whether the temperature dependence

96

A. Vernov, W. Steele / Simulation of N 2 adsorbed on graphite I

i

!

I

i

0 45 90 135

4~ 90I 135180 '

4,

4, •

I

!

l

I

I

'

V

I

~1'

I

"~1

0 45 90 135

I

I

I

45 90135

4,

I

I

I

q,

!

4,

!

I

i

i

I

90 135180

|

I

I

I

45 9 0 135 0 45 90 135

4,

I

45

4, i

04590135

I

045 90155

'~5

4,

45 90 135180

4,

Fig. 10. Plots of the distribution of in-plane orientation angle ~ are shown for a number of simulations. The commensurate lattice at 25 K (# 1) has lost its herringbone ordering in both the first and second layers (left- and right-hand parts of the figure, respectively). The X-compressed systems shown in the second row (:~4, 5, 6) are ordered, but show a significant loss in order in going from 25 to 35 K. Distributions for the Y-compressed case are shown in the bottom row for 25 K (~8) and 35 K (#9).

of the potential energies given there is close to that for a classical h a r m o n i c 1 oscillator of ~ T . per degree of freedom. This a m o u n t s to 5 / 3 2 = 0.16 for the change from 25 to 35 K. The changes in m o l e c u l e - s o l i d energy a m o u n t to 0.47 a n d 0.21 ( c o m m e n s u r a t e ) ; 0.54 a n d 0.24 (X-compressed) a n d 0.63 a n d 0.25 (Y-compressed) for the first- a n d second-layer molecules, respectively. (The n u m b e r s for the X- a n d Y-compressed layers should be the same, if the proposed r e a r r a n g e m e n t of Y- into X- is correct. A p r o b l e m in this context is that it is likely that the r e a r r a n g e m e n t was not completed d u r i n g the equilibration period a n d consequently, calculated averages for Y-compressed systems include configurations that are not yet in the X-compressed form.) T h e modes of m o t i o n that can affect the molecule solid interaction are primarily v i b r a t i o n a n d libration in the p l a n e p e r p e n d i c u l a r to the surface; however, the d e p e n d e n c e of this energy on position in the plane parallel to the surface is non-negligible for first-layer molecules and thus may well affect the t e m p e r a t u r e d e p e n d e n t averages. Indeed, two degrees of freedom yield a value

A. Vernou, IV. Steele / Simulation of N 2 adsorbed on graphite l

I

I

97

I

#I

u

-I.O

J o

u.

i

I

I

I

0

-l.O cos ~ l

I

I0 I

#4

-I.0

-I.(

J

0

I

I

I

0

-I.0 cos /3

i

-I .0 COS ,~ I

I

I 0

I

I

I

G

1.0

1.0

Fig. 11. Distributions for the out-of-plane orientation angle /~ at 25 K are shown for first-layer (left-hand panel) and second-layer (right-hand panel) molecules. Dashed lines indicate the distribution for monolayers and solid curves are for bilayers. A commensurate (~1), and two X-compressed systems, # 4 (small compression) and # 7 (large compression).

of 0.32 for the change from 25 to 35 K, which only amounts to 2 / 3 of the first-layer value, so one concludes that the periodic terms in the energy are significant but do not contribute the full harmonic oscillator value to this " h e a t capacity". The small values of = 0.23 obtained for the change in second-layer energies may be understood by realizing that these molecules are no longer near the bottom of the molecule-solid potential well so that averages over their harmonic motions do not give rise to average energies of ½ kT, but something much smaller. (In the limit of a linear variation in molecule-solid energy near the second-layer distances, the average energy is constant at all temperatures.) The temperature dependence of the N 2 - N 2 in-plane energy amounts to 0.27 and 0.20 (commensurate); 0.27 and 0.15 (X-compressed); and 0.28 and 0.27 (Y-compressed) for the first and second

98

A. Vernov, IV. Steele / Simulation of N 2 adsorbed on graphite

layer, respectively. In all cases, these values are much lower than the value for two in-plane vibrational plus one in-plane librational degree of freedom. The basic problem here is that these layers cannot be considered to be strictly two-dimensional; that is, N z - N 2 potential energies are obtained by averaging over configurations whose probability is governed in part by the molecule-solid and the interlayer NR-N 2 interactions. Consequently, the simple picture leading to ½kT of thermal excitation per degree of oscillatory freedom is no longer valid for the intralayer potential energies. Now turning to the energies obtained for T = 25 K, one sees that even the molecule-solid Um*~,~varies from one simulation to another. For example, the addition of a second layer causes of decrease in U,*,,~ equal to 0.29 (commensurate) and 1.05 (large X-compression). This is due to the change toward coplanarity in the out-of-plane orientation of the first-layer molecules when a second layer is added. Since the molecule-solid potential is quite sensitive to this orientation near its minimum, the observed energy changes are not unexpected. The relative stability of the layers simulated here depends upon their free energies. Nevertheless, the total potential energies can give some indication of which systems might be stable. (Of course, the transformation of Y-compressed into X-compressed is a proof, but transformation involving density changes or interlayer transfers, such as compressed ~ commensurate, do not occur on the timescale of the simulations.) At 25 K, the commensurate monolayer has an energy lower than the (greatly) X-compressed layer by 0.65 reduced units (runs Z~3 and 6), but this difference is much smaller for the two bilayers (runs # 1 and 7) where AU~*~= 0.11. It is straightforward to perform a minimum energy calculation for these layers (in effect, at 0 K). Assuming that all molecules are coplanar with the surface and that the incommensurate phases are perfectly regular (no long-wavelength periodic variations in density), the energies of the monolayers were calculated as a function of the lattice spacing for equilateral triangular, X-compressed films with Y = Y(commensurate) and Y-compressed films with X = X(commensurate). Since the molecule-solid potential is unchanging in this calculation (except for a completely commensurate film), it is not included. Curves of the N 2 N2 energy for these monolayers are shown in fig. 12. It is evident that the Y-compressed layer is always the least stable; that the commensurate layer would be unstable relative to either an X- or to a uniformly-compressed (equilateral) layer, except for the added energy due to the periodic part of the molecule-solid potential. For these systems, this periodic energy amounts to - 0 . 4 3 . Its inclusion would give the commensurate layer an energy very nearly the same as the minimum energy X-compressed layer, which has a 7% linear dimension change. At 25 K, our results for the commensurate (run :~ 1) and the X-compressed with a 6% dimension change (run ~ 6 ) confirm the picture of a low-energy uniaxial layer, since U~] is lower by 0.3 (0.6 at 0 K) relative to

A. Vernoo, W. Steele / Simulation of N 2 adsorbed on graphite

-I

\

I

99

I

Y

e

-13

-

1.2

I

I

I. I

I.O

0.9

P / Pcomrn. Fig. 12. The N 2 - N 2 potential energy for monolayers of coplanar molecules is shown as a function of density (i.e., of lattice spacing) for triangular lattices which are: equilateral (equi.); uniaxial, with spacing along one cartesian axis held at the commensurate value but a variable spacing along either the X-axis ( X ) or the Y-axis. See fig. 8, bottom for definition of the axis orientations and an example of a Y-compressed layer. In all cases, the in-plane orientation angle was varied to find the minimum energies shown.

commensurate. However, U~,] is higher by 1.46, due partly to the periodic potential (0.43) and partly to orientational effects. (The use of periodic boundary conditions will suppress any long-wavelength density variations.) The evidence produced here supports the experimental finding that a first-order commensurate ---, X-compressed transition occurs for N 2 monolayers; however, the experimental X-dimension appears to be = 3% compressed at ---20 K, which is significantly smaller than the calculated 7% at 0 K. Nevertheless, the calculated energies are quite close to one another and thus, a careful calculation of the entropic contribution to the free energies is needed before a definitive result can be obtained for the thermodynamics of the model monolayers. It is interesting that the second-layer N2-N 2 energies tend to be more negative than for the underlying first layer (excluding the Y-compressed systems). Since the densities of the two layers have been fixed to be equal (which is not necessarily the case at equilibrium), the most obvious difference between the layers is that molecules are more nearly coplanar in the underlayer. It would appear that the strongly negative energies available to N 2 N 2 pairs in non-parallel configurations give rise to the decrease in second-layer energy relative to first. In the coplanar herringbone structure found for the low-temperature monolayers, the angle q5 made by the molecular axis with the lattice glideline

100

A. Vernov, W. Steele / Simulation of N 2 adsorbed on graphite

60 °

I

I

50" ,p 40 c

7"~3// / 30 ¢

6-.t. 0

1.2

I

I

I.I

I.O

0.9

P/ Pcomm Fig. 13. The in-plane angles for the herringbone lattice of N 2 configurations corresponding to the minimum energies of fig. 12 are shown by the solid lines. The values of ~b given are those between the axes and the crystal glidelines, so that the two sublattices differ only in the sign of ~p. Also shown are most probable values of ~ obtained from plots of probability distributions for 25 K such as those in fig. 10; the run numbers for each point shown are indicated in the figure.

is a variable quantity which depends upon temperature, compression and, in simulations, upon the details of the N 2 - N 2 potential energy function. For 0 K, these angles are generated as part of the energy minimization routine. Fig. 13 shows the variation in the minimum energy value of in-plane angle with compression; also shown are the most probable values obtained at 25 K from plots such as those in fig. 10. Values for the commensurate systems are omitted because they are orientationally disordered at this temperature. The calculation show that increasing X-compression and increasing temperature have the effect of making the most probable orientation of the molecular axes be more nearly parallel to the glidelines. The interpretation of this behavior is greatly complicated by the simultaneous changes in out-of-plane orientation. 4. Conclusions

It is perhaps not surprising to find that molecular motions in bilayers at these low temperatures appear to be nearly harmonic. Nevertheless, the

A. Vernov, W. Steele / Simulation of N 2 adsorbed on graphite

101

changes in mean orientation angle observed and in addition, the shifts in molecular position that occur in the Y-compressed ~ X-compressed transition indicate clearly that the motion in these systems is not as restricted as implied by a harmonic oscillator description. We conclude that orientational and translational equilibrium is reached within each layer, even though the limited time duration of the simulations makes it highly unlikely that the interlayer transfers required to achieve overall equilibrium will occur. In spite of this limitation, it is notable that the two layers are strongly coupled in the sense that imposition of a second layer causes the first-layer molecular axes to undergo a significant shift toward coplanarity with the surface. Furthermore, the second layers appear to all be "in registry" with the first-layer lattices and, when orientational ordering is present, have preferred axes for orientation which are coupled to those for the first-layer molecules. All these coupling effects are present even though the amplitudes of librational and translational oscillatory motion are significant. Of course, the second-layer molecules feel the molecule-surface potential less strongly and consequently, their orientational ordering relative to the surface is less pronounced than for the first (either in monolayer or bilayer systems). When the out-of-plane orientational ordering distributions show a significant degree of non-coplanarity, it should be emphasized that the distributions have a tendency toward being bimodal. This can be interpreted as a mixture of (nearly) coplanar and perpendicular molecules, rather than a gradual change of the most probable value of the angle fl as the entire layer reorients. The model showing such a mixture is the pinwheel, where perpendicular N 2 molecules are at the centers of hexagonal arrays of flat molecules. Such pinwheels can occasionally be seen visually in snapshots of the instantaneous configurations, but attempts to quantitatively characterize their presence were unsuccessful. Finally, our results are not helpful in answering questions concerning the wetting properties of N 2 multilayer films at T < 40 K. As noted above, the transfer of molecules from the initially postulated structures to some other, be it bulk crystal or even layers with different densities than initial, appears to be too slow for the attainment of equilibrium in these simulations. (The transfer rate does increase with increasing temperature and is readily observed at 75 K [21]; however, the simulations give fluid films at that temperature with complete wetting behavior, but at least partial roughening.) Even a direct evaluation of free energy in the simulation would not be possible unless and until the layer densities are at their equilibrium values. On the other hand, the finding that the orientations and positions of the second-layer molecules are strongly correlated with those in the underlayer is helpful in understanding the experimental work on nitrogen films at these temperatures. For example, Zhang, Kim and Chan [28] have found that the partly completed bilayer film consists of monolayer plus two-layer parts which behave as coexisting phases

102

A. Vernov, W. Steele / S i m u l a t i o n o f N 2 adsorbed on graphite

with definite thermodynamic properties. Our results are certainly consistent with these data and indeed, the experiments suggest that simulations at even higher coverages migh help in understanding layering and wetting behavior at these temperatures. Acknowledgements Support from grant DMR-8419261 from the the National Science Foundation is gratefully done while AV was an exchange visitor at Penn International Research and Exchange Board Education of the USSR.

Material Research Division of acknowledged. This work was State under the auspices of the and the Ministry of Higher

References [1] Y. Larher, J. Chem. Phys. 68 (1978) 2257. [2] T.T. C h u n g and J.G. Dash, Surface Sci. 66 (1977) 559. [3] S.C. Fain, Jr., J.F. Toney and R.D. Diehl, in: Proc. 9th Intern. Vacuum Congr. and 5th Intern. Conf. on Solid Surfaces, Madrid, 1983, Ed. J.L. de Segovia (Imprenta Moderna, Madrid, 1983 pp. 129-137. [4] J.K. Kjems, L. Passell, H. Taub, J.G Dash and H.D. Novaco, Phys. Rev. B13 (1976) 1446. [5] J. Eckert, W.D. Ellenson, J.B. Hastings and L. Passell, Phys. Rev. Letters 43 (1979) 1329. [6] R.D. DieM, J.F. Toney and S.C. Fain, Jr., Phys. Rev. Letters 48 (1982) 177. [7] R.D. Diehl and S.C. Fain, Jr., J. Chem. Phys. 77 (1982) 5065. [8] R.D. Diehl and S.C. Fain, Jr., Phys. Rev. B26 (1982) 4785. [9] J. Piper, J.A. Morrison, C. Peters and Y. Ozaki, J. Chem. Soc. Faraday Trans. 79 (1983) 2863. [10] M.H.W. Chan, A.D. Migone, K.D. Miner, Jr. and Z.R. Li, Phys. Rev. B30 (1984) 2681. [11] K.D. Miner, Jr., M.H.W. Chan and A.D. Migone, Phys. Rev. Letters 51 (1983) 1965. [12] Q.M. Zhang, H.K. Kim and M.H.W. Chan, Phys. Rev. B32 (1985) 1820. [13] C. Peters and M.L. Klein, Mol. Phys. 54 (1985) 895; Phys. Rev. B, submitted. [14] J. Talbot, D.J. Tildesley and W.A. Steele, Mol. Phys. 51 (1984) 1331. [15] J.P. Joshi and D.J. Tildesley, Mol. Phys. 55 (1985) 999. [16] J. Talbot, D.J. Tildesley and W.A. Steele, Chem. Soc. (London) Faraday Disc., in press. [17] A.D. Migone, H.K. Kim, M.H.W. Chart, J. Talbot, D.J. Tildesley and W.A. Steele, Phys. Rev. getters 51 (1983) 192. [18] M. Nielsen, K. Kjar and J. Bohr, J. Electron. Spectrosc. Related Phenomena 30 (1983) 111. [19] K. Morishige, C. Mowforth and R.K. Thomas, Surface Sci. 151 (1985) 289. [20] J. Talbot, D.J. Tildesley and W.A. Steele, Surface Sci. 169 (1986) 71. [21] A.V. Vernov and W.A. Steele, Langmuir, in press. [22] W.A. Steele, Surface Sci. 36 (1973) 317. [23] W.A. Steele, J. Phys. (Paris) C4 (1977) 61. [24] W.E. Carlos and M.W. Cole, Surface Sci. 91 (1980) 339; G. Vidali and M.W. Cole, Phys. Rev. B29 (1984) 6736. [25] C.S. Murthy, K. Singer, M.L. Klein and I.R. McDonald, Mol. Phys. 41 (1980) 1387. [26] D.J. Evans and G.P. Morriss, Computer Phys. Rept. 1 (1984) 297; Phys. Letters 98A (1983) 433. [27] R.M. Lynden-Bell, J. Talbot, D.J. Tildesley and W.A. Steele, Mol. Phys. 54 (1985) 183. [28] Q.M. Zhang, H.K. Kim and M.H.W. Chan, Phys. Rev. B32 (1985) 1820.