A molecular dynamics simulation of tje uniaxial phase of N2 adsorbed on graphite

Surface Science 169 (1986) 71-90 North-Holland, Amsterdam

A MOLECULAR DYNAMICS SIMULATION OF N 2 ADSORBED ON GRAPHITE

71

OF THE UNIAXIAL PHASE

J. T A L B O T a n d D.J. T I L D E S L E Y Department of Chemistry, The University, Southampton S09 5NH, UK and W.A. STEELE Department of Chemistry, The Pennsylvania State University, University Park, Penn.~vlvania 16802, USA Received 2 August 1985; accepted for publication 7 November 1985

A computer simulation study of a nitrogen monolayer adsorbed on a graphite substrate at low temperatures is reported. The adsorbed phase was slightly compressed relative to the commensurate 7'3 x v/3 phase. The compression was taken to be uniaxial; i.e. a 5% change in the large intermolecular spacing along a glide line of the oriented herringbone structure (UXI phase). Thermodynamic properties were evaluated together with 0rientational and translational ordering parameters at two temperatures, one above ( = 40 K) and one below ( = 15 K) the in-plane disordering transition. For purposes of comparison, simulations of the commensurate phase at these temperatures are reported, together with a 15 K simulation of a phase that has been uniaxially compressed in the direction perpendicular to that of the UXI phase (UYI phase). The simulations indicate that the UXI phase is stable but the UYI phase tends to transform into domains of UX1; it is also concluded that the compression necessary to form the UXI phase from the commensurate does not produce a significant change in the out-of-plane ordering, at least at the lower temperature, but does bring about changes in the in-plane ordering of these molecules.

1. Introduction N i t r o g e n f o r m s a v a r i e t y of t w o - d i m e n s i o n a l p h a s e s w h e n a d s o r b e d o n the basal p l a n e of g r a p h i t e . A t t e m p e r a t u r e s b e l o w 27 K a n d c o v e r a g e s less t h a n o n e m o n o l a y e r , the n i t r o g e n f o r m s an o r i e n t a t i o n a l l y o r d e r e d , ~ 3 × v/3 c o m m e n s u r a t e solid, C, w h i c h coexists with a d i l u t e t w o - d i m e n s i o n a l gas [1-4]. T h i s solid u n d e r g o e s a r o t a t i o n a l p h a s e t r a n s i t i o n at 27 K [5,6] b u t r e m a i n s c o m m e n s u r a t e up to the m e l t i n g t e m p e r a t u r e , w h i c h is d e n s i t y d e p e n d e n t [7]. If c o v e r a g e is i n c r e a s e d a b o v e o n e c o m m e n s u r a t e m o n o l a y e r the d i f f r a c t i o n d a t a i n d i c a t e that the N 2 m o l e c u l e s f o r m a t r i a n g u l a r i n c o m m e n s u r a t e solid ( T I ) w i t h an i n t e r m o l e c u l a r s p a c i n g of 0.404 n m in c o n t r a s t to 0.426 n m in the C phase. T h e c e n t r e of mass s t r u c t u r e is the s a m e as that f o u n d for a d s o r b e d k r y p t o n at high d e n s i t i e s [8]. A t low t e m p e r a t u r e s the T I solid is likely to be 0 0 3 9 - 6 0 2 8 / 8 6 / $ 0 3 . 5 0 © Elsevier Science P u b l i s h e r s B.V. ( N o r t h - H o l l a n d Physics P u b l i s h i n g D i v i s i o n )

J. Talbot et al. // Molecular dvnamic.~ w m u l a t . m

72

orientationally ordered, with a 2-out herringbone (or "pin-wheel") structure [9a]. However, Fain and Diehl [gb] have also discovered a new uniaxial incommensurate phase (UI) which exists over small range of coverage for densities between that of the C and TI solids. Although diffraction has onl>. been observed at one temperature (15 K) this phase is likely to exist over a range of temperatures, and below approximately 27 K it is orientationally ordered. The existence of this phase has been confirmed by synchrotron X-ray diffraction [10], and the orientational disordering of the UI solid has been observed by calorimetry [11]. The structure determined by LEED is shown in fig. lb, and may be contrasted with the corresponding (7 structure shown in fig. la. Both the UI and C structures consist of two sublattices of molecules arranged in a herringbone pattern. The UI phase is obtained by compressing the C phase in one direction towards the glide-line, marked A in fig. 1. The extent of the linear compression is estimated from LEED t o range tip to a maximum of 5%, while the lattice spacing in the direction orthogonal to the compression remains constant at the value of the C solid. The U I structure can exist in one of three separate domains, with the glide-line A pointing in the direction shown in fig. lb or at an angle of + 120 ° to this direction. In reality, these domains coexist on the graphite surface. There remain a number of questions to be asked about the new phase. The glide-line is a line of translation-reflection symmetry. Its presence in the Ul phase (inferred by systematic spot absences in the LEED pattern) means that the angle between the crystal axis and the N 2 bond is the same for molecules in both sublattices. The LEED experiment does not determine the value of this angle. Furthermore the experiments tell us little about the distribution of out-of-plane molecular orientation and its change with changing surface density. Finally it is not clear why a compression should occur uniquely in a direction towards the glide-line A. At first glance, an orthogonal compression A t

0 426nrn I

0 426nmi

0 369nm ia)

0 352 nm (b)

Fig. 1. (a) T h e ~/3 x V- c o m n l e n s u r a t e s t r u c t u r e for N• a d s o r b e d o n g r a p h i t e . T h e naolcculen a r c o r i e n t a t i o n a l l y o r d e r e d in the l w o - s u b l a n i c e h e r r i n g b o n e s t r u c t u r e . (b) T h e u n i a x i a l i n c o m m e n s u r a t e s t r u c t u r e in w h i c h tile o r i e n t a t i o n a l s t r u c t u r e is r e t a i n e d a n d the m o l e c u l e s arc c o m p r e s s e d by 5% in the d i r e c t i o n p e r p e n d i c u h i r to line A. A is a glide-line. (There is ~l s e c o n d glide-line p e r p e n d i c u l a r to A not s h o w n in the figure.)

J. Talbot et al. / Molecular dynamics simulation

73

where molecules in the same sub-lattice are more closely packed appears to be a reasonable alternative UI structure. In a recent paper, we demonstrated that the C phase can be simulated using the molecular dynamics technique [12]. By performing a series of simulations at constant coverage we located the transition between the orientationally ordered and disordered phase C phases for our model of the N 2 layer. The agreement between simulation and experiment was satisfactory and in this paper we have used a similar model to examine the UI phase. In particular we will adopt the center of mass (COM) structure determined by the LEED experiment and will study the in-plane and out-of-plane orientational ordering. In section 2 we describe the potential model used in this work. The N2-N 2 interaction has been slightly modified from that used in our first study [12] by adopting an effective pairwise model which has been tested against a wide range of solid and liquid properties. The simulation technique has been described in detail elsewhere [12] and in section 2 we comment on the additional problem of simulating an incommensurate solid. In sections 3 and 4 we describe the translational and orientational structure of the low and high temperature UI phases. The paper offers a comparison between the properties of the UI and C phases at the same temperature, and to this end we have performed two commensurate phase simulations using the new model potential. In section 5 we discuss a simulation study of the compression orthogonal to that observed in the LEED experiment, and section 6 contains our conclusions. Peters and Klein [13,14] have also studied the compressed N 2 overlayer on graphite by Monte Carlo simulation. These workers assume somewhat different intermolecular potentials than those used here; in addition, they studied the two-phase (or striped domain) region for a large number (2704) of molecules in one case and the one-phase uniaxially compressed system [13] to a higher compression (9%) than that chosen in this work (5%). These simulation studies are generally consistent with those to be presented here; detailed comparisons will be given below.

2. The potential model and simulation The N 2 N 2 interaction used in this work is the X1 model of Murthy et al. [15]. This is a two-center Lennard-Jones model with an additional quadrupole quadrupole interaction. In this simulation the quadrupole ON, is represented by placing partial charges along the molecular axis. (A charge of q = 2 0 N , / I ~ , at each nitrogen atom, and - 2 q at the center of the bond.) These partial charges give an accurate representation of the quadrupole, but fail to model the higher electrostatic moments accurately [17]. There is still some uncertainty in the actual value of the N 2 quadrupole moment. Accurate ab-initio calculations at the CI level give a value of -5.51 × 10 4o C m 2 [18]

74

J. Talbot et aL / Molecular dynamics simulation

whereas birefringence experiments give - 4 . 6 x 10 4o C m 2 [19]. In our previous simulation of N~ adsorption we used a value of - 5 . 0 7 x 10 4o C m 2. The X1 model uses a value of - 3 . 9 1 × 10 4l) C m 2, wlaich represents a 15c/ reduction from the experimental value. The X1 model q u a d r u p o l e should be thought of as an effective q u a d r u p o l e designed to include other elements of the overall charge d i s t r i b u t i o n such as exchange forces and higher electrostatic moments. The X1 model is preferred to the earlier two-center L e n n a r d - J o n e s models [20] as it gives a good fit to the properties of both bulk solid and liquid N,. In particular it gives an accurate lattice energy for the solid and reproduces the translational p h o n o n frequencies, although it consistently over-estimates the frequencies of the librational modes in the c~-solid. The low value of 0~<, in the X1 model may be useful in modelling adsorption since the use of this potential in the simulation will reduce the transition temperature from the 33 K found in our previous study [12] towards the experimental value of 27 K [5], However we should be cautious in adjusting a potential model to reproduce an experimental transition temperature because of the difficulties of d e t e r m i n i n g this temperature accurately for the model, given the small system size in the s i m u l a t i o n a n d the complication of degenerate d o m a i n structures. We have used the F o u r i e r - e x p a n d e d molecule surface potential due to Steele [21]. Two terms in the expansion are retained and c o n v e n i e n t expressions are given in the appendix of ref. [12]. The molecule surface potential parameters are changed slightly from those used in ref. [12] since they are d e t e r m i n e d from the N : - N 2 potential by mixing rules. The isostcric enthalpy at zero coverage for the present model is 9.98 kJ reel ~ which lies between two recent experimental d e t e r m i n a t i o n s (9.7 kJ reel 1 [22] and 10.4 kJ nlo] i [23]). The parameters used in the molecule molecule and moleculc surface potential are summarized in table 1. A similar potential has recently been used by Bruch [24] in a study of N 2 adsorption. In this paper we have not included substrate mediation of either the dispersion or the electrostatic interaction between N~ molecules. "laDle 1 Paranlctcr~, used in the modelling of N~ adsorbed on graphite. The NN paralllewrs arc model X1 of ref. [13]. The ('(" parameters are suggested by %teele[14] and the NC parameters arc determined by the Lorent~,. Berthelot mixing rules ~..,/k u oN~ ON, IN, %< /'k)~ a( (

(K) (nm) (104o (" m~) Into)

36.4 0.331 3.91 {).1()9~

(K)

2~.0

(nm)

(N(/k.

(K)

o N~

(nm)

0.340 ~1.92 0.336

J. Talbot et al. / Molecular dynamics simulation

75

The simulations were performed on 98 N 2 molecules in a three-dimensional cell. The cell thus has dimensions of 4.92 x 2.98 nm and is periodic in the x and y directions, it also has a reflecting wall which is 2.0 nm above the surface. The equations of motion are solved using a 4th order predictor-corrector method and the orientations are handled using quaternions. Because of the relatively large size of the two-dimensional cell, forces and torques could be accurately and rapidly calculated using the minimum image convention without spherical cutoff. The energy is conserved to better than one part in 105 over the complete run, using a timestep of 2.5 x 10 -~5 s. Runs described in this paper consisted of 2 x 104 timesteps after an equilibration of 4 x 103 timesteps. Configurations were stored on tape at intervals of 10 timesteps. In a molecular dynamics simulation of an adsorbed layer, the box has to be chosen so as to contain an integral number of unit cells of the substrate. Thus a molecule which crosses the boundary of the central cell and enters as its image from the opposite side experiences no discontinuity in surface forces. The structure of the adsorbate in an incommensurate solid is unlikely to fit into the same simulation cell as the surface lattice unless there is an accidental degeneracy of the structures at specific distances. Fortunately the UI solid, at 5% compression, has a repeat distance in the x-direction perpendicular to the glide-line A such that this choice of cell dimension allows periodic boundary conditions for both the substrate and adsorbate. In the UI phase this MD cell contains an even number of rows thus accommodating a possible two sublattice structure. In fig. 2 we show a starting configuration for the UI solid. The molecules are arranged at angle of 45 ° to the glide-line A in this initial configuration. We have performed two simulations of the UI phase and two of the C phase at

Y

",,

',,"

.\.

.

:\i

[ i',i

F i g . 2. T h e two

.',

k

"

"

;,

i

" <

"

" \"

'

'\'

.

~

.

.

~.

.

.

.

. \.

.

.',,.

i

~" i

: "~ ~ i',i

:

i',i

.

i 5, i

initial configuration

sublattice

initially

"

structure•

at an angle

The

"~

for the simulation glide-line

points

o f ~ r / 4 t o t h i s line• T h e

of the UI solid. 98 molecules

along spacing

the space-fixed between

rows

y-axis

and

are organized the molecules

in t h e x - d i r e c t i o n

is 0 • 3 5 2

in a are nm.

J. Talbot et al. / Molecular dynamws simulation

76

Table 2 Average energies evaluated for the four simulations of the (' and UI phases using the model discussed in section 2 (see table 1). The standard error in the temperature is given in parentheses. The total energy E is conserved and there is an error of +0.02 kJ tool l associated with the configurational energies. U~ is molecule-surface contribution and Um is the molecule molecule contribution. Note that E = U, + U,,, + K , where K is the average kinetic energy ( 2.5RT. for N2)

UI ("

UI (,

T K)

E (kJmol i)

U, (kJmol i)

I4,' (kJmol i)

14.75 +0.21) 16.18 -+0.23) 40.25 ( + 0.35) 40.19 ( + 0.38)

- 12.23

9.68

2.85

12.13

9.75

2.71

10.89

9.15

2.57

- 11/.96

9.32

2.47

a p p r o x i m a t e l y the s a m e t e m p e r a t u r e . T h e final t e m p e r a t u r e s of the s i m u l a t i o n s are s h o w n in table 2, t o g e t h e r with the c o n f i g u r a t i o n a l e n e r g i e s c a l c u l a t e d o v e r the run. As we expect, the c o m p r e s s i o n of the p h a s e c a u s e s the m o l e c u l e surface e n e r g y to b e c o m e less a t t r a c t i v e as m o l e c u l e s m o v e out of registry. T h e c o m p r e s s i o n also c a u s e s an i n c r e a s e in the a t t r a c t i o n b e t w e e n N , m o l e c u l e s . In the low t e m p e r a t u r e solid the c h a n g e in N~ N~ a t t r a c t i o n is d o m i n a n t and the U I solid has a slightly l o w e r c o n f i g u r a t i o n a l e n e r g y t h a n the (" solid [ - 1 2 . 5 3 kJ mol i ( U I ) , - 1 2 . 4 6 kJ tool i ((7)]. In the high t e m p e r a t u r e solids the c h a n g e in m o l e c u l e s u r f a c e i n t e r a c t i o n d o m i n a t e s , a n d o n e o b t a i n s [ 1 1.72 kJ tool i ( U i ) . 11.79kJmol i(C)].

3. Translational structure T h e t r a n s l a t i o n a l s t r u c t u r e of the a d s o r h a t e can be d e t e r m i n e d by m o n i t o r i n g two o r d e r p a r a m e t e r s , if ()PI gives the c o r r e l a t i o n b e t w e e n the p o s i t i o n s of the a d s o r b a t e m o l e c u l e s and the s t r u c t u r e of tile u n d e r l y i n g s u r f a c e hittite, we h a v e OP~ = ~

,<

~,

I

cosk~.r,+cosk ]

2.r,

)

.

1)

If OP~ d e t e r m i n e s the c o r r e l a t i o n b e t w e e n n e i g h b o r i n g a d s o r b a t e m o l e c u l e s . ~sc can write OPt=

N 21 2

<

Y', ,v Y', c o s k ~ - r , ; + c o s k , - r , ~ 1t 1 "

/)X.

{2)

J. Talbot et al. / Molecular dynamics simulation

77

For the C solid, k 1 and k 2 a r e the reciprocal lattice vectors for the vf3 × c o m m e n s u r a t e lattice. Thus, k 1 = (27r/a)(2~3, 0) and k 2 = (27r/a)( - 1 / 3 , l / v / 3 ) , where a = 0.246 nm is the unit of distance used in the simulation. In eq. (2), ~ and r,j are two-dimensional vectors representing the projection of positions and separations onto the surface. In the case of the UI solid the reciprocal lattice vectors are k 1 = (2~r/a)(0.7, 0) and k 2 = ( 2 ~ r / a ) ( - 0 . 3 5 , 1 / ¢ 3 ) . O P 2 is well-defined, but OP 1 is simply the correlation between the position of the centers of mass of the N 2 molecules and their initial positions at the beginning of the simulation. The values obtained for OP l and OP 2 are given in table 3. It can be seen that OP 2 is high for all four simulations; indicating that as the temperature is increased to 40 K both solid structures are preserved. (The slight fall in order parameter with temperature is characteristic of increased oscillation around the lattice sites.) If the solid remains fixed with respect to the underlying surface lattice then OP~ -- (OP 2)1/2.

(3)

This relationship describes the C solid accurately and indicates that in this phase the molecules are strongly associated with the m i n i m u m energy adsorption site at the center of the graphite hexagon. However in the case of the UI solid OP 1 is frequently smaller than (OP2) 1/2 and tends to vary during the simulation. This variation in OP~ is caused by the floating of the UI solid over the graphite. The center of mass of the adsorbate moves in a direction perpendicular to the glide-line. As one row of molecules leaves the low energy sites in the center of the graphite hexagons, another row of atoms becomes " c o m m e n s u r a t e " and there is only a small energy barrier to lateral movement. The direction of the movement must be influenced by the initial assignment of Table 3 The translational order parameters for the adsorbed solid. The numbers in brackets are the standard error calculated in sub-averages formed from 500 timesteps. For the uniaxial solid OP 1 is not constant during the run and the values given in the table are averages calculated over 500 timesteps at 1 × 104 and 2 x 104 timesteps, z is the most probable height for COM of the N 2 molecule and Zm,,x is the maximum height observed in the simulation T (K)

OP I

OP 2

UI

14.7

0.090 ] -0.332/~

0.907 (_+0.002)

C

16.2

0.889 ( -+0.002)

UI

40.2

0.940 ( -+ 0.001 ) 0.845 ] 0.592 /

C

40.1

0.856 ( _+0.004)

0.739 ( -+0.008)

0.732

Z (nm)

Zmax (nm)

0.332

0.381 +_0.012

0.332

0.381 +0.012

0.332

0.455 -+0.012

0.332

0.4305 + 0.012

J. Talbot et al. / Molecular dvnamtcs simulatton

78

10 ~ ° 2, ;, o

m

O

O

Q

Q

D

O

O

4

Q

O

O

aYO.O 0

-10

1JO

15 K B/10 3

2JO

2'5

Fig. 3. ()PI as a f u n c t i o n of t i m e d u r i n g the s i m u l a t i o n . T h e values are s u b - a v e r a g e s c a l c u l a t e d over the p r e c e d i n g 500 timesteps. (0) UI solid. 14.7 K. (zx) U1 solid, 40.2 K. K B = n u m b e r of (imesteps.

the velocities to the molecules. In the case of the low temperature U I solid we have observed floating along the x-axis, both in a positive and in a negative direction. In fig. 3 the values of OP~ calculated over groups of 500 timesteps are plotted. The order parameter initially falls rapidly from a value 0.94 and becomes negative: at 20 × 103 timesteps it appears to be rising again towards zero. The projections of the C O M of the N 2 molecules from this simulation are shown in fig. 4a. At timestep zero row 1 is along the edge of the box in the c o m m e n s u r a t e position. In fig. 4a, we show only the projections for the last 104 timesteps of the run. The total lateral movement is indicated by the arrows above the cell. By the end of t h e simulation row 1 is again commensurate. Because of the periodic b o u n d a r y conditions in the simulation, row 13 may eventually cross the edge of the cell, at which point O P 1 will have regained a value of 0.94. The motion shown in fig. 4a is highly cooperative, with an average distance between rows which is essentially constant at 0.351 nm throughout the simulation. (This contrasts with a spacing of 0.369 nm in the C solid.) Interestingly, this p h e n o m e n o n appears to be less p r o n o u n c e d at high temperature. The plot of OP~ for a high temperature UI simulation is shown in fig. 3. There is a region of lateral movement between 17 × 10 ~ and 24 × 10 ~ timesteps but the solid returns to its original position with respect to the underlying surface by timestep 25 × 103, In other simulations of the high temperature UI phase, we have observed a slow fall in OP 1 but this occurs after 104 timesteps without any movement in the adsorbate. A possible explanation is that the m o v e m e n t in the high temperature U1 phase is likely to be less

J. Talbot et aL / Molecular dynamics simulation

~J'

a

~: ~..~i "

: ~i "

~/~

~." i i'~! i ;~'. ~ i~i ] i~ : ~'~ ~ i*-i ~ ~"i i ?'~: .

i~ "

-t~-

i ~':

~

• ak~-"

' ,,b~ "

~

"

• , ~

79

,~

¢o.,' -

'

"

'

"

~

"

~ i ~. i ~ i

"

~

,~.

' ,~:,

~

"

"

"

"~,

-

~-"

'

"

'~4

.

, _..

" ~ .:_.

'

'~,'

'

"~

"

i i~'. ~ ?,,~ ~ .,,,

'~

ok,~

"

"

',~-'

"

'.~'

"

~*"

"

",~'

'

",~

",~'

%~-'

, ~ ,

"

'

,

- ~

(1)

•

(13)

b i • . . i~i . . . . i . . ~--i

/

!~i

]

~:

~

i

: ~ "

~i .~/.

i ~- • i q". . . . . . . . . ia~'

: ~'~'i

~:

i ~.

: ~.

. ,~.

~ ~ ~ . + _ _

~:

~,,:

-'~

.,~

~-

~ ~

~,~ .~

~:

i

'~ .~ ~

~

i

/

Fig. 4. T h e t r a j e c t o r i e s o f the C O M o f t h e N 2 p r o j e c t e d o n t o t h e s u r f a c e . T h e p r o j e c t i o n s s h o w the m o t i o n f r o m t i m e s t e p 7 × 103 to t i m e s t e p 2 × 104. (a) U ] solid, 14.7 K. T h e a r r o w s m a r k t h e l a t e r a l m o v e m e n t o f c o l u m n 1 o v e r the e n t i r e s i m u l a t i o n . (b) U1 solid, 40.2 K.

cooperative since the molecules are already vibrating substantially around their equilibrium position and that the row-by-row commensuration associated with the low temperature movement cannot be achieved. The projections from a high temperature UI run are shown in fig. 4b. (Note that these simulations are performed in the presence of a static external field. Although the energy is conserved, there is no absolute requirement for the conservation of linear momentum. During the course of the equilibration net m o m e n t u m is removed at each timestep but this condition is relaxed during the run to ensure that projections are generated in the microcanonical ensemble. Since the external field is symmetric we do not expect net displacement in the center of mass on a long timescale.) Turning now to structure in the dimension perpendicular to the surface, the most probable height for finding a molecule above the surface, i.e. the

80

J. Talbot et a L / Molecular dynamics simulation

m a x i m u m in n ( z ) , is found to be constant at 0.3321 _+ 0.012 nm for the temperatures considered .here. The m a x i m u m height at which the center of mass of a molecule is found, increases with temperature. It is also slightly higher for the UI solid than C solid at 40 K. This effect is small and there is no evidence for second layer promotion or significant out-of-plane structural rearrangement caused by the compression.

4. Orientationai ordering 4.1. In-plane In-plane orientational ordering is described using the distribution functions n ' ( 0 ) and n"(q,). These distribution functions apply to the two sublattices, and are proportional to the n u m b e r of molecules with in-plane orientations between q, and g, + dq,. ~ is the angle between the projection of the intermolecular axis onto the surface and the glide-line A. The herringbone ordering can be established by monitoring the order parameters

I
OP4 =N-S

E sin 20, i=1

f +"/2sin

2q,

n'( q, ) dO

~,~/2 r,/2

J

(4)

n'(dt) ) dO ~r/2

/ + ' ¢ / 2 s i n 2¢) n " ( ~ ) dO k

N,~

sin 2~5; = ;=1

1.7/2 ~rr/2

where the sums are over the N ' and N " molecules in sublattices one and two. If all the molecules are aligned at 45 ° to the glide-line, OP 4 = o p s - 1 : if the distribution is uniform, the order parameters fall to zero. (This is a more natural choice than the OP 4 and OP 5 defined in ref. [12], where the angles are measured with respect to the side of the simulation cell in an alternative ground-state structure of the herringbone.) The value of OP4 and o p s obtained in this work are shown in table 4. As expected, in both the UI and C solids the order parameters from both sublattices are equal to within the error in the simulation. In the high temperature phase OP4 and o p s are zero to within the estimated error, indicating substantial in-plane disordering. Examination of the n ( ~ ) distribution functions shows limited residual ordering (of a different symmetry than the UI phase) which is not monitored by the order parameters. OP 4 and ops. Since a more complete description of orientational ordering is contained in the distribution functions, we define

n(o) =

[,,'(,) + n"(,)],

(6)

81

J. Talbot et al. / Molecular dynamics simulation

Table 4 The orientational order parameters [eqs. (4), (5), (10)]; the numbers in brackets are the standard error calculated from sub-averages formed from 500 timesteps

T

oP,

oP5

oP~

0.925 ( _+0.001) 0.810 ( _+0.002) - 0.018 ( _+0.04) 0.027 ( _.+0.02)

- 0.924 ( _+0.001) - 0.807 ( + 0.002) 0.010 ( _+0.03) - 0.013 ( _+0.02)

-

(K) U1 C U! C

14.7 16.2 40.2 40.1

0.496 ( _+0.001 ) 0.470 ( _+0.001) 0.261 ( +_0.002) 0.335 ( _+0.002)

which has the useful p r o p e r t y that

n ( ~ ) = n(-q~).

(7)

The function n(q~) can be written as a F o u r i e r series n(q~) = ½ a 0 + ~

a; cos 21q~.

(8)

/=l

B e c a u s e n(q~) is a s y m m e t r i c a b o u t ~r/2, o n l y c o s i n e t e r m s of a r g u m e n t 2l o c c u r in the series. T h e sin t e r m s d i s a p p e a r b e c a u s e o f the e q u i v a l e n c e o f the t w o s u b l a t t i c e s (we h a v e c h e c k e d that this is the case in the s i m u l a t i o n p r e s e n t e d here). Values o f the a; c o e f f i c i e n t s for the f o u r s i m u l a t i o n s are s h o w n in table 5 a n d the n(q3) f u n c t i o n s in fig. 5. T h e l o w - t e m p e r a t u r e solids clearly s h o w two s u b l a t t i c e structures. In the C solid n(q,) is r e a s o n a b l y s y m m e t r i c a l a r o u n d ~ r / 4 as e v i d e n c e d by the high v a l u e of a z and the low v a l u e of a 1 (a~ s h o u l d be z e r o for a f u n c t i o n w h i c h is s y m m e t r i c a r o u n d ~r/4). In the U I solid, the m a x i m u m in n ( q 0 is shifted t o w a r d s s m a l l e r angles a n d the a~ t e r m is an o r d e r of m a g n i t u d e h i g h e r t h a n in the C solid. T h e h i g h e r F o u r i e r c o e f f i c i e n t s of the U I solid are larger, i n d i c a t i n g that the series is m o r e slowly c o n v e r g e n t , c o n s i s t e n t w i t h the m o r e s h a r p l y p e a k e d n(,~). Thus, the u n i a x i a l c o m p r e s s i o n

Table 5 The coefficients in the Fourier expansion of n(~) T

a 0

a I

a 2

1.998 2.000 2.000 1.999

0.256 0.026 0.200 0.001

- 1.468 1.091 -0.069 -0.005

a 3

(/4

(/5

0.579 0.239 -0.018 0.012

0.368 0.080 0.006 -0.009

(K) UI C UI C

14.7 16.2 40.2 40.1

0.516 0.133 0.159 0.119

82

J. Talbot et al. / Molecular dvnamw,~ ,~irnulation

UI 14-7K

r

i

4-

C

u c

16.2K

2A 0

,

o

~/2

¢ 4

Ul 4 0 2 K

b

0 4 [

401K

0

½

Fig. 5. The in-phme orientational distribution function n((/,) showi~ for thc four simulallow, gixcn in table "~

has sharpened the in-plane distribution while also pushing the molecular o r i e n t a t i o n s toward the glide-line. The average value of ~ is

1

,/"2

n ' ( O ) dO

fT-//F( O dO

f ~"20 n " ( O ) dO t9)

f_........~~- 2n"( ~ ) dO

J. Talbot et al. / Molecular dynamics simulation

83

can be easily estimated from the simulation. It is 42.80 ° + 0.15 ° for the C solid and 41.0 ° + 0.15 ° for the uniaxial solid at 17 K. (This value of 42.80 ° differs from the value of 4 4 . 5 ° + 1.5 ° given in ref. [12], a result of the difference in the potential functions used in the simulations, in particular the lower value of 0N2 used in this work. Unfortunately there is no experimental estimate of this angle available.) In the high temperature C phase n(q,) has six-fold symmetry (over the range 0-27r), as evidenced by the high value of a 3 relative to the other coefficients. This weak angular ordering in the temperature phase is due to the crystal field of the adsorbate. Interestingly, in the UI solid, the six-fold symmetry is less pronounced than in the C phase (note the high value of a~). From the graph of n(q,) for the high temperature UI phase, it is clear that the six-fold symmetry is lost by the suppression of the peaks at _+~r/2. In the compressed solid the molecules cannot point as easily in the direction of the compression. 4.2. O u t - o f - p l a n e distributions

The out-of-plane ordering is described using the distribution n(cos fi). This describes the probability of finding a molecule with a tilt of fl away from the normal to the surface. The degree of tilt can be determined by monitoring the order parameter

op~=~

1

P2(cos fl,)

/=1

(10)

-

/

n(COS/~) d c o s f l

op~ = - 0 . 5 if the molecules are all parallel to the surface and 0 if they are randomly oriented. The values of OP 3 are shown in table 4. The negative number obtained in all simulations indicate that the most probable orientation is parallel to the surface. The compression to form a uniaxial layer has little effect on OP 3 in the low temperature solids. The n(cosfl) distributions obtained for the low temperature UI and C solids are almost identical and similar to that shown in fig. 12b of ref. [12] so that we do not reproduce them here. It is important to realize that the molecules alter their in-plane ordering rather than tilting away from the surface, for the compression used. In the high temperature solids the compression in the UI solid causes more tilting towards the normal than in the corresponding commensurate solid. (Compare rows three and four in table 4.) The distributions n(cos/3) for the high temperature solids are shown in fig. 6 where it can be seen that n(cos fl) for the UI solid has a more pronounced up-turn around +~r/2 and - ~ ' / 2 than the distribution for the C solid. These features in n(cos/3) are due to the transitory appearance of pin-wheel-like configurations during the run. The pin-wheel is the favored structure for an ordered two-dimensional triangular

84

J. Talbot et al. / Molecular dynamics simulation

32•

UI 40.2 K

u

C 40.1K

2

o

O0

0!5

cos(I])

1!0

Fig. 6. Distributions of cos/3, where/3 is the angle between the molecular axis and the normal to the surface plane, are shown for both high temperature phases. (Symmetry requires that these functions be the same for positive and negative values of cos ,8.)

lattice of quadrupoles in the absence of a strong field at the crystal surface [27]. One molecule at the center of the pin-wheel is oriented perpendicular to the surface while the six neighboring molecules are parallel to the surface and

. . . . . .

.

• I .

.

.

.

~-

.

: ":

t

~

.

.

.

'

'I

.'

.

.

. .'.

.

r

.

.

.\

.

.,"

.

:

.

.

.

.

L/:

, - : -

.

~,

.

:

'

.'~

"

.\

,t'.....,./

.

.

\

-

•

. i .

.

.

•

.

.

.

.

•

.

/.

.

.

/:

•

,

.

,'~

'

."

~"

".

'\

'

•

"

.

.

.

:1:

:

:\:

:

(

:

.

.

.

:

.

,

.

./,

,

.

,

.

."

I

'

:~ :

.

.

•

•

:/:

.

.

.

./

.

:

: ,!

a

.

:

:

.

.

.

:l :

jt

'

.

.

-

~ '

Fig. 7. A snapshot of a molecular configuration in the high temperature uniaxial phase. Projections of the molecular axes on the surface lattice are shown, and the circle indicates a molecule that is upright. Note that the six neighbors to this molecule form an excellent representation of a pin-wheel structure• (Other such pin-wheels can be located elsewhere in the snapshot, but none with a vertical molecule in the center.)

J. Talbot et aL / Molecular dynamics simulation

85

arranged so that their molecular axes are at a fixed angle relative to each other. We do not expect perfect pin-wheels to form in the UI lattice since it is not equilateral triangular, but transient local pin-wheels are clearly visible. A configuration of the UI high temperature simulation is shown in fig. 7 and a pin-wheel center is marked with a dashed circle. It seems likely that these structures dominate in the triangular incommensurate phases which occur as the compression is increased. We have not examined the temperature dependence of the rotational phase transition in detail for the UI solid. This would require a substantial number of simulations at a constant coverage of 6.72 nm -2. This transition certainly occurs in the model since we are able to observe ordered and disordered UI phases. However, a simulation of 98 molecules would be too small to make an accurate appraisal of the order of the transition.

5. Compression parallel to the glide-line As adsorbed N 2 undergoes a uniaxial contraction towards the glide-line, the angle that the molecules make with respect to the glide-line decreases and the molecule-molecule energy becomes more attractive. It is of interest to examine what happens during a compression parallel to the glide. In the C solid, the average spacing between the centers of N 2 molecules in the same sublattice, i.e. along the direction of the glide-line, is 0.426 nm. In this simulation this distance is reduced to 0.3977 nm, which gives a 6.7% increase of density. Use of a cell of 2.95 n m × 5.96 nm containing 120 molecules allows for periodic replication of the adsorbate and the underlying lattice with the possibility of a two sublattice herringbone structure. The initial configuration is shown in fig. 8. (It is difficult to construct a cell which is more nearly square without doubling the number of N 2 molecules in the simulation.) A simulation of this U YI phase was performed using 4000 steps of equilibration followed by 104 steps to calculate the equilibrium averages. (The proposed phase is named according to the direction of the compression.) The final temperature of the run is 17.55 + 0.55 K and the total configurational energy is - 12.35 kJ mol-1. The structure is slightly less energetically favorable than the UI solid ( U = - 1 2 . 5 3 kJ mo1-1) although much of this difference may be due to the small differences in temperature between the two simulations. Since the simulated heat capacity of the commensurate phase is -- 45 J mo1-1 K l [12], correction to the energy for A T = 2 K would be = 0.09 kJ mol 1, giving a residual difference which is significantly less than the uncertainty in either simulated energy value. The reciprocal lattice vectors of the UYI solid are 2~/a(-1/3, 15/14v/3) and 2~r/a(2/3, 0). At T = 17.55 K, the calculated values of OP1 and OP 2 are 0.56 and 0.60 respectively. These low values indicate a significant structural

86

J, Talbot et al. / Molecular dynamics simulation " :,~:" :" :~: " :" :;,: " :" :~,: •

. / ,

•

. / ,

.

. / .

•

.~.....~...-.~,.-...;,. ...>.....;...) ... .,:,... - . k . ' . . . & . -...<. • .;'.-.

. z . .

• . . . y . . .

•

,

-

/

•

.

.

•

f

..

.

.

4

.

.

,':'...

"

7

"

"

2

'

,

' f

'

•

"

' / "

•

. " .

'J"

• ".

• ~,'

. • .

.'~"

, .

• • .

'

"#'"

.

.

",t'"

."

.

' ( ' .

. ' 4 " . " .

,/.

.

,/.

•

,j.

: :,,:

• j.

-

7

"

'

/ "

-

"

•

%.¢ '

• "

.

. •

.~,

"('

. • .

...,~.

•

:,,:

-

.

.:,.....,~....:,....~..

Y

•

.

\'

.,,.....~,......~, .

.

•

/

"

"K'

.

.

,

, 4 . . .

'

. " .

'

,

'

• , , ,

. / . . .

•

•

.2..

" ~ ' . ' .

.

.,-.

•

-j.

: :x:

.

: :'.,: .

-'~X

Fig. 8. A pictorial representation of the initial configuration chosen for the UYI crystal; i.e., for the nitrogen layer compressed in a direction parallel to the y axis shown on the figure.

I

a 4 q~

Z

2

I

-I.0

-018

-0.6

-0.4

-0.2

0,0

0.2

0.4

016

0.8

cos(B) i

i

i

b z

C' -

80

- 40

0

40

80

Fig. 9. Orientational distributions obtained for the U Y ] phase at = ]5 K. O u t - o f - p l a n e orientations are given by n(cos f l ) and the in-plane distribution by n ( ~ ) .

J. Talbot et aL / Molecular dynamics simulation

87

rearrangement from the starting configuration. Furthermore the value obtained for OP 3 is -0.409, which is significantly lower than the value of - 0 . 4 7 0 found for the UI and C solids at low temperature. The plot of n(cos r ) shown in fig. 9a shows more pronounced out-of-plane angular ordering with quite distinct peaks at + 1. The fact that these peaks are uneven indicates that the one end of the molecule remains pointing away from the surface and is indicative of long-lived and stable pin-wheel structures. It is also found that OP 3 and OP 4 have fallen to low values of 0.215 and - 0 . 2 5 0 respectively, indicating loss of in-plane orientational order• The apparent loss of in-plane ordering is confirmed by the n ( ~ ) distribution in fig. 9b. This unexpected result can be explained with reference to fig. 10 which shows the last configuration from the simulation. This snapshot indicates that the y compression has caused a structural rearrangement. The new structure consists of two domains of herringbone which are incommensurate. The clear pin-wheel structures which have been circled in the diagram form a natural grain boundary between domains of herringbone. There is a mismatch of structure at the top and bottom edge of the box resulting in a small region of orientational disorder. It can now be seen that the parameters OP 3 and OP 4 are small because the herringbone ordering is rotated relative to the initial glide-line along the space-fixed y-axis which is used as origin in the calculation of OP 3 . . . ',-". '~. : : • .5.

~:. . . < ... ' -' ;•,-. '-- .

• ..:..;.-...~.

• .<.....~. ....~....~.

..~.

.'~.'.'-I"•

•

"-.~.

•

.

~

"..~.. " .I

\.

~: - :

~

.

~.

"

'

~ " " "/'"

..~.

..r."

"

,~.

." . \ . '

./,..

~'

-:

:

~

,

•

:

~i

.~.

.

• " " ~:'

..~.

" - ~) - "

" ". ' ~ "

"I ; - , ' '

I "'"

~. r . "

,.. • . • .~ , , - ' ~ ~ ~ "I . .: . . . . I. ~

::: . . : ~ . ~

.

.

• I",'.~.

. k.. .--.

.

".~..

.' ~.

• ' . / - ' I.

"~" -

•.i.

..~...

".'.'

• ' •

:.:

"

-. ~- .,..,. ~ ) . . ]. .

.~.

:

.

• " " " I: . :". t""~'.t "

~' " - 1 ' " I'

l - ' -' ' ' ; - '' ' ' ' ' I' ' "

-~ ."

I

".#..

" ~:

. \.

..l.....

' .;--

" :~" ~ " • ~ • " :'~'-"

' ."

-.,;.

•.~.

-;."

~...

" ~:"' '~ k~2 " 7 ~:" " " '/ " ~-

..(.

..~.

•

-

. ~ . . . I .

. . . . ~ . . . - . : . . . . ~ . .

p ~

~

'.,:.'-.f•

~. ..~...._ " .,;-.

"~ : ~ :

~~'~'~' I .'."~ ~

..1

..~....

: ~

~

~

,~ " - .<. '-. -"~: -"- - 4 - - .

" . ' - ~ ' " . " J - " i -

;

"&'"

"

' .'~.

" ,"

.r

."

•

~-- ,,~. • i- ,~. --. ,~. • I "I" .'"T' "I'.:~

: '-. - ~ . ' .

.I

Fig. 10. A snapshot of a molecular configuration reached in the simulation of the UYI phase. On comparing with figs• 2 and 8, it can be seen that the molecules have reoriented to form two of the three UXI sublattices, separated by a disordered region containing several local pin-wheel structures. (For clarity, this figure shows the computation cell plus an image.)

88

J. Talbot et al. / Molecular dvnamics simulation

and OP 4. (This is the only simulation in which we have observed the transformation from one herringbone structure to that belonging to another domain.)

6. Conclusions Much remains to be done if one is to give a detailed characterization of these uniaxially compressed N 2 monolayers by computer simulation. For example, the experiments suggest that the transition from C to UI is first order, at least at low temperature, with the usual finite density and energy differences between the coexisting phases. Clearly, many more simulations at different densities and temperatures would be needed to establish this point. In view of the difficulty in obtaining proper matching at the periodic boundaries, simulations over the requisite range of lattice parameter may not be feasible for reasonable numbers of molecules. Furthermore, these boundary conditions have a tendency to favor the ordered structures that match the translational periodicity and may thus yield metastable phases in the relevant density region. Nevertheless, several interesting conclusions emerge from the simulations reported here. First, the UYI phase is clearly unstable relative to the UX1 system. This instability is relieved by orientational changes from UYI to a mixture of two of the UXI sublattices (rather than by translational rearrangements see above). Of course, it would be helpful to have estimates of the free energies, or more precisely, of the free energy differences between these dense monolayer phases. However, these differences are quite small, as indicated by the size of the energy differences obtained in this study. In comparing these results with experiment, one must take account of an additional complication, which is that the magnitude of the periodic terms in the model gas-solid potential is not very well known. Both theory [25] and simulation [26] indicate that the actual barriers to translational across the surface may be as much as two times larger than those used here and indeed, Peters and Klein [13,14] have used a molecule-solid potential with a larger periodic component than that used here. In the absence of direct comparison, it is hard to estimate the effects of an increase in the periodic terms, but it is worth noting that a considerable alteration in the two-dimensional melting point of N: on graphite could be produced by changing these terms [26]. Clearly, any commensurate ~ incommensurate transition will be quite sensitive to this part of the potential, and major changes in the predicted phase equilibria cannot be excluded when the energy balance is as delicate as that observed in this study. As for as the N 2 N 2 potentials go, these are constrained by the requirement that they give realistic values for the properties of the bulk material. From that point of view, the differences between our choice and that of Peters and Klein are not large. On the other hand, the value of N= quadrupole moment (which is

J. Talbot et al. / Molecular dynamics simulation

89

proportional to the partial charge on an N atom) used here is significantly smaller than that employed in our earlier study of the commensurate N 2 layer. In fact, in unpublished work, we have observed that a decrease in partial charge of 20% causes the orientational transition temperature for the commensurate layer to drop from ---33 to ---23 K. Thus, the orientational properties are quite sensitive to the electrostatic part of the potential. For example, only minor adjustments are needed to achieve agreement with the experimental or the Peters-Klein transition temperatures of 27 K. The effect of such an adjustment on the uniaxial-commensurate-triangular-incommensurate phase equilibrium is hard to guess. It is unlikely that any of the known direct simulation methods of evaluating chemical potential will be sufficiently accurate to give useful results for the dense low temperature systems of interest here. Approximate alternatives such as assuming harmonic oscillator solids and evaluating free energy from the calculable phonon frequency distributions are also dubious because of anharmonicity, particularly for the in-plane N 2 reorientations. An additional problem that arises in considerations of phase equilibria in these systems is whether any of the cases studied by computer simulation is an equilibrium phase in the sense of having the lowest free energy consistent with the areal density and temperature. For example we have shown that a y-compressed layer is unstable but we have not located the most stable phase, which could be x-compressed, uniformly compressed or even compressed in both x and y directions but by different amounts. Since the results will depend upon temperature and density as well as the details of the N 2 - N 2 and Nz-surface potential and quite possibly, the system size and shape, definitive answers would appear to require an enormous computational effort. Nevertheless, these simulations indicate that it costs very little energy to form the UXI phase from the lower density commensurate phase. This uniaxial phase does exhibit an order-disorder transition in the in-plane orientational degree of freedom which appears to be similar to the corresponding change in the commensurate phase, in agreement with experiment. The fact that a 5% linear compression can be achieved in the simulation without greatly altering the degree of co-planarity of the N 2 molecules and the surface is a result that cannot be easily obtained by experiments on the real system. Finally, the spontaneous appearance of pin-wheel structures at or near what appear to be crystal domain walls is an unexpected way of relieving the strains associated with such defect structures.

Acknowledgments JT acknowledges a SERC studentship, and DJT and WAS a travel grant from N A T O (216.80). Support by a grant from the Division of Materials Research of the N S F is also acknowledged.

90

J. Talbot et al. / Molecular dynamics simulation

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

J.K. Kjems, L. Passell, D. Taub, T.G. Dash and A.D. Novaco, Phys. Rev. BI3 (1976) 1446. J. Eckart, W.D. Ellenson, J.B. Hastings and L. Passel[, Phys. Rev. Letters 43 (1979) 1329. R.D. Diehl, M.F. Toney and S.C. Fain, Jr., Phys. Rev. Letters 48 (1982) 177. R.D. Diehl and S.C. Fain, Jr., Surface Sci. 125 (1983) 116. A.D. Migone, H.K. Kim, M.H.W. Chan, J. Talbot, D.J. Tildesley and W.A. Steele, Phys. Rex'. Letters 51 (1983) 192. N.S. Sullivan and J.M. Vassieu, Phys. Rev. Letters 51 (1983) 658. Y. Larher, J. Chem. Phys. 68 (1978) 2257: M.H.W. Chan, A.D. Migone, K.D. Miner and Z.R. Li, Phys. Rev. B30 (19841 2681. S.C. Fain, Jr., M.D. Chinn and R.D. Diehl, Phys. Rev. B21 (1980) 4170. (a) H. You and S.C. Fain, Faraday Disc. 80 (July 1985); (b) R.D. Diehl and S.C. Fain, Jr., Phys. Rev. B26 (1982) 4785. M. Nielsen, K. Kjaer and J. Bohr. J. Electron Spectrosc. Related Phenomena 30 (1983) 111. Q.M. Zhang, H.K. Kim and M.H.W. Chan, Phys. Rev. B32 (1985) 182(/. J. Talbot D.J. Tildesley and W.A. Steele, Mol. Phys. 51 (1984) 1331. C. Peters and M.L. Klein, Mol. Phys. 54 (1985) 895. C. Peters and M.L. Klein, Phys. Rev., in press. C.S. Muthy, K. Singer and I.R. McDonald, Mol. Phys. 41 (1980) 1387. W.A. Steele, The Interaction of Gases with Solid Surfaces (Pergamon, New York, 1974) sects. 2.8, 4.9. C.S. Murthy, S.F. O'Shea and I.R. McDonald, Mol. Phys. 50 (1983) 531. R.D. Amos, Mol. Phys. 39 (1980) 1. A.D. Buckingham, R.L. Disch and D.A. Dunmur, J. Am. ('hem. Sot. 90 (1968) 3104: N.J. Bridge and A.D. Buckingham, J. Chem. Phys. 40 (1964) 2733. P.S.Y. Cheung and J.G. Powles, Mol. Phys. 32 (1976) 1383; 30 (1975) 921. W.A. Steele, J. Phys. (Paris) C4 (1977) 61. Y. Grillet, F. Roquerol and J. Roquerol, J. Phys. (Paris) C4 (1977) 57. J. Piper, J.A. Morrison, C. Peters and Y. Ozaki, J. Chem. Soc. Faraday Trans. 1, 79 (1983) 2863. L. Bruch, J. Chem. Phys. 79 (1983) 3148. G. Vidali and M. Cole, Phys. Rev. B29 (1984) 6736. Y.P. Joshi and D.J. Tildesley~ Mol. Phys. 55 (1985) 999. A.B. Harris, O.G. Mouritsen and A.J. Berlinsky, ('an. J. Phys. 62 (1984) 915.