A Monte Carlo simulation study of Nitrogen on LiF(0 0 1)

Applied Surface Science 256 (2010) 2974–2978

Contents lists available at ScienceDirect

Applied Surface Science journal homepage: www.elsevier.com/locate/apsusc

A Monte Carlo simulation study of Nitrogen on LiF(0 0 1) A.K. Sallabi a, J.N. Dawoud b,*, D.B. Jack c a

Department of Physics, 7th October University, P.O. Box 2666, Misurata, Libya Department of Chemistry, The Hashemite University, Zarqa 13115, Jordan c Department of Chemistry, University of British Columbia-Okanagan, Kelowna, BC V1V 1V7, Canada b

A R T I C L E I N F O

A B S T R A C T

Article history: Received 8 October 2009 Received in revised form 19 November 2009 Accepted 19 November 2009 Available online 2 December 2009

The adsorption of N2 gas on the LiF(0 0 1) surface is studied by canonical Monte Carlo (CMC) computer simulation. These pffiffiffiresults pffiffiffi show that N2 forms an ordered structure where the molecules are arranged in a unit cell of pð2 2  2Þ R 45 symmetry at temperatures below 23 K with 50% coverage. The nitrogen molecules are tilted by 538 from the surface normal and have the same azimuthal orientation along diagonals, with diagonals alternating their orientation. Beyond 23 K, the molecules become azimuthally disordered but with residual short-range order. No change in the position of the peak of the polar (tilt) angle distribution was observed above the transition temperature. This transition is purely of the order– disorder type. ß 2009 Elsevier B.V. All rights reserved.

Keywords: Physical adsorption Monte Carlo simulations Phase transition LiF

1. Introduction In the last twenty years, experimental techniques have been improved and made a wealth of information on adsorbate structures at low temperatures, Epitaxial growth and 2D phase transitions, exhibited by molecules adsorbed on solid surfaces [1– 13], available. This was attracting considerable attention for both experimentalists and theoreticians. One of the interesting systems that widely investigated is the N2/NaCl [2] system. Previous Monte Carlo simulations of the adsorbed molecules systems were able to reproduce most of the experimental results [1–3], in addition the simulations of the N2/NaCl [2] show that this system undergo an order–disorder phase transition with nonuniversal critical exponents and fall in the same universality class as the XY model with cubic anisotropy. The values of the critical exponents depend on the strengths of an anisotropic external potential provided by the substrate (NaCl). Note here order/ disorder phase transition was also observed for CO/NaCl(0 0 1) system using classical Molecular Dynamics simulations [14]. In order to extensively test the modern theory of phase transition [6], a different system with different anisotropic symmetry from N2/NaCl is needed. The best candidate is the N2/ LiF system. The same conditions as that of N2/NaCl case have been applied to simulate the N2/LiF system. However, because of its smaller lattice constant, N2 molecules are expected to be adsorbed

* Corresponding author. E-mail address: [email protected] (J.N. Dawoud). 0169-4332/$ – see front matter ß 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.apsusc.2009.11.060

at every other Li+ site instead of at every Na+ site as in the case of the NaCl substrate. In this paper we report on a set of the Metropolis Monte Carlo (canonical ensemble) simulations of the N2/LiF(0 0 1) system to investigate possible stable structures, molecular orientations and phase transitions. 2. Potentials and methods The methods and potentials used in this work are similar to those applied in the simulations of phase transitions in the CO2/ MgO [1], N2/NaCl [2] and CO/MgO [3] systems. Therefore, an outline of the interaction potential of N2/LiF(0 0 1) system will present here. The interaction potential that applied here is divided into two parts. The first part describes the interaction between the nitrogen molecules adsorbed on the surface, in which, the intermolecular potential is partitioned into contributions, V mm ¼ V elec: þ V re p: þ V dis p:

(1)

where, the three terms are electrostatic, repulsion, and dispersion contributions, respectively. In our model, the nitrogen molecule is represented as two atoms that interact individually with the atoms of other N2 molecules on the surface as well as with the ions of the surface below. Therefore, the intermolecular potential was calculated, within two-body approximation, as a pair-wise sum of atom–atom interaction potential. The electrostatic interaction was created by assigning point dipoles of equal magnitude but opposite sign, located on each atomic sites of N2 molecule, and both pointing towards the center of the bond. Therefore, the

A.K. Sallabi et al. / Applied Surface Science 256 (2010) 2974–2978 Table 1 Summary the values of the various experimental properties of N2 multipoles in addition to their molecular polarizabilities that required for electrostatic model of the N2 molecule.

generated at the surface by the substrate. As applied elsewhere [1], the calculation of the N2-surface electrostatic interactions was performed using the advantage of the periodicity of the substrate ** (lattice symmetry) and hence, the electric field E ðr i Þ above the (0 0 1) face of a fcc ionic crystal can be determined from the * electrostatic potential Fðr Þ of the surface as follows,

N2 Dipole (D) Quadrupole (D, A˚) Ocutpole (D, A˚2) Hexadecapole(D, A˚3) Polarizability (a1) (A˚3) Polarizability øø (A˚3) Polarizability ? (A˚3) (a2) (A˚5)

0.0 1.41 [15] 0.0 3.01 [16] 1.74 [17] 2.19 [18] 1.513 [18] 77.96 [16]

**

5 X Ci j

2n

where, *

Fðr Þ ¼ 

ms Velec: ¼

where, the parameters Ai j and bi j , are the Born–Mayer parameters that characterize the strength and inverse range of the repulsion interaction between ith and jth atoms located on different N2 molecules, r i j is the distance between nitrogen atomic sites ij distributed on different molecules, and C2n represents the London dispersion coefficients. The values of the Born–Mayer parameters for N atom were calculated using the following formulas [21,22],   ai þa j 1 bi þb j ; and bi j ¼ Ai j ¼ ðbi þ b j Þ  e bi þ b j

V

C8 (A˚8kcal/mol)

C10 (A˚10kcal/mol)

N2

1010 [19]

9354 [20]

–

 2  X * * *2 1 1 ==*2 E ðr i Þ  m i  ai? E ? ðr i Þ  ai E == ðr i Þ 2 2 i¼1

TT

ðr i j Þ ¼ Ai j  expðbi j  r i j Þ 

5 X n¼3

(4)

f 2n ðr i j Þ

ij C2n

! (5)

ðr i j Þ2n

ij where, Ai j and bi j are the Born–Mayer parameters, and the C2n dispersion coefficients represent the strength of mutually induced attractive forces due to electron correlations. The f 2n ðr i j Þ is the damping function that is given by;

" f 2n ðr i j Þ ¼ 1 

2n ðb  r Þk X ij ij k¼0

k!

#  expðbi j  r i j Þ:

(6)

According to Tang–Toennies potential model, it is necessary to estimate the repulsion parameters and the long-range dispersion coefficients that characterize the interaction of the N atom with the ions of the surface. These parameters can be constructed using combining rules from data for N–N and like–like ion interactions. Below an outline of the procedure that used to calculate these parameters is presented. The like–like ion repulsion parameters, Li+–Li+ and F–F, are taken from references [4] and [5]. These values along with the parameters of the nitrogen atom–atom interaction, presented previously, are used to construct the different atom–ion Born– Mayer parameters using the Gilbert and Smith combining rules [21,22]. In particular, the Born–Mayer parameters of two different

Table 2 Summary the values of the dispersion coefficients of molecular N2. C6 (A˚6kcal/mol)

"   #    

z exp  2p 4e 2px 2py pa ffiffiffi þ cos cos a 1 þ expð 2  pÞ a a

where, the sum is over the two atomic sites of nitrogen molecule and the induction energy is proportional to the atomic polariz== abilities perpendicular ðai? Þ and parallel ðai Þ to the N2 molecular axis. Tang–Toennies potential [25] with summation over two-body interactions has been applied to describe the repulsion and dispersion terms of molecule–surface interaction using the following formula,

where, ai ; a j ; bi and b j are the nitrogen atomic softness parameters that are estimated using the WMIN program [23]. The dispersion coefficients of atomic nitrogen have been estimated from the dispersion coefficients of molecular nitrogen that are presented in Table 2 using mathematical combining rules [24]. The values of the repulsion–dispersion parameters are calculated and listed in Table 3. The second part of the potential is used to describe the total interaction between the adsorbed N2 molecule and the ions of the LiF substrate. This kind of potential is considered to be the sum of classical electrostatic interactions, plus contributions arising from repulsion and dispersion interactions. They are assumed to be pairwise additive so that the molecule–surface potential may be written as the sum of two-body interactions between each constituent atom of the admolecule and the individual ions of the substrate. Using the electrostatic model of a nitrogen molecule * described above, the two point dipoles m i , located on atomic N * sites at positions r i will interact with the electric field E ðr i Þ

Molecule

(3)

where, a = 2.85 A˚ is the like-ion separation on LiF lattice, and e is the magnitude of the electronic charge on an individual ion. Furthermore, this electric field will also interact with dipole moments induced on the atomic sites of the molecules by the electric field itself. The electrostatic plus induction interaction energy between a nitrogen molecule and the substrate may be calculated using the following formula,

(2)

r 2n n¼3 i j

*

E ðr Þ ¼ rFðr Þ;

electrostatic energy may be estimated as the sum of dipole–dipole interactions. The values of the atomic dipoles are chosen so as to reproduce the various experimental molecular properties of N2 multipoles [15–18] that are summarized in Table 1. The repulsion–dispersion of molecule–molecule interaction was calculated using Buckingham potential formula, Virej p:dis p: ðr i j Þ ¼ Ai j  eðbi j ri j Þ 

2975

Table 3 The calculated repulsion–dispersion parameters of the N2–LiF and N2–N2 interaction. Interaction

Aij (kcal/mol)

bij (A˚1)

C6 (A˚6kcal/mol)

C8 (A˚8kcal/mol)

C10 (A˚10kcal/mol)

N(atom)–Li+ N(atom)–F N–N

45325.4 23915.2 31000

4.58 3.32 3.48

13.055 256.63 252.8

51.4 1555.635 2340

247.905 11551.31 49540

A.K. Sallabi et al. / Applied Surface Science 256 (2010) 2974–2978

2976

species may be calculated using the following formulas,

bi j ¼

2  bii  b j j bii þ b j j

c

and

Ai j ¼

ðAii  bii Þ  ðA j j  b j j Þ

bi j

d

(7)

where c ¼ b j j =ðbii þ b j j Þ and d ¼ bii =ðbii þ b j j Þ. The resulting values of Ai j , and bi j for a N atom interacting with the ions of the surface are listed in Table 3. The values of the dispersion coefficients, C 6 ; C 8 and C 10 , for a N atom interacting with the surface ions, Li+ and F, have been calculated using combining rules that derived from the Casimir– Polder theory of the dispersion coefficients [26]. For example, the C 6 coefficient was calculated using the following relation, C6i j ¼

2  C6ii  C6j j  ai  a j C6ii ða j Þ2 þ C6j j ðai Þ2

discussed conveniently with respect to this coordinate system by defining in the usual way a polar angle u (tilt with respect to the surface normal) and an azimuthal angle w (angle between the x axis and the projection of the molecular axis onto the surface plane). During the MC simulations an ensemble of 72 adsorbed molecules are placed in the surface potential of the ionic substrate, where the size of the central square box was 34.2 A˚ (12 surface square lattices). Periodic boundary conditions in the lateral directions (x and y) were imposed as well as a cutoff radius of 13.5 A˚ for the molecule–molecule interactions. Although no restrictions were placed on the ‘‘motion’’ of the molecules in the z direction (away from the surface) the molecules nonetheless stayed in the vicinity of the surface due to the attraction by the surface potential and the cold temperatures at which the simulations were run.

(8) 3. Simulation results

where, i refers to a N atom, j refers to either a Li+ or F ion, and a is the dipolar polarizability of the corresponding atom or ion. The C 8 coefficients are calculated using a lengthy procedure that is explained in detail elsewhere [25], whereas the C 10 coefficients of the N (atom)–Li+(F) interactions have been estimated using the following relation [27],

In agreement with experiments [12], MC simulations at 1 K confirm that a single molecule physisorbed on LiF(0 0 1) surface

2

ij C10 

49  ðC8i j Þ

40  ðC6i j Þ

(9)

The values of the atom–ion dispersion coefficients that applied in this work are listed in Table 3. We have used the Metropolis Monte Carlo simulation (canonical ensemble or Q(N,V,T)) since it allows for the sampling of a large number of possible configurations at nonzero temperature. The number of molecules N, the volume V and temperature T are fixed during any simulation. In our simulation, the surface is not allowed to move, whereas, the adsorbed molecules are considered to be rigid (fixed bond length) and subjected to a random translation in all directions (lateral and perpendicular) as well as rotations around the nitrogen atom. The amplitudes of the moves were adjusted independently and maintained a 50% acceptance probability. For each Monte Carlo cycle a randomly chosen molecule is moved to a new random position or orientation. The energy of the new configuration is calculated and the difference in energy with the old one, DE, is evaluated. The Metropolis sampling algorithm is applied to decide whether to accept or reject the new configuration. Application of this procedure yields configurations that are weighted according to the Boltzmann distribution and hence are suitable for studying average structure, energies, and other computed properties of the system. The results of the MC simulations are presented in the next section. In experiments related to our work, molecules are deposited under UHV conditions onto the (1 0 0) face of cleaved single crystals of ionic surfaces. Thus it is reasonable to assume that large sections of these surfaces are defect free and hence can be represented as an infinite periodic surface. We assume that the substrate can be regarded as a semi-infinite solid with the surface atoms occupying the positions they would in a bulk system. This allows for a coordinate system to be affixed to the crystal surface. The origin of the coordination system is chosen to be located on an anion site in the plane of the surface with x and y axes running along [1,1,0] and [1,1,0] crystallographic directions respectively. The z axis runs along the [0,0,1] crystallographic direction normal to the surface. It is worth noting that for reasons of simplifying the simulations, the x and y axes are rotated by 458 around the surface normal with respect to the principal axes of the bulk. By taking advantage of the periodicity of the ions of LiF(0 0 1) surface, the periodic surface potential may be represented as a two dimensional Fourier series. The orientation of the N2 molecules can be

Fig. 1. The final configurations of a monolayer of N2 on LiF(0 0 1) surface at 1 K (top panel), 20 K p (middle ffiffiffi pffiffiffipanel) and 25 K (bottom panel). The monolayer at 1 K forms an ordered pð2 2  2Þ R 45 structure. The solid white circle represents the Li+ site, whereas the brown ones represent the F sites.

A.K. Sallabi et al. / Applied Surface Science 256 (2010) 2974–2978

with a binding energy of 1.07 kcal/mol (experimental value is 1.2  0.1 kcal/mol)[12]. Our results show that N2 molecule prefers to adsorb perpendicular with respect to the surface plane on top of a Li+ cationic site at a height of 2.825 A˚. At a 50% coverage, where 72 molecules of N2 are placed on a (12  12) patch of LiF surface that contains a 144 Li+ sites and an equivalent number of F sites, the Monte Carlo simulations show the existence of an ordered phase structure at 1 K. Inpthis ffiffiffi phase, pffiffiffi the N2 molecules arranged in a unit cell of pð2 2  2Þ R 45 symmetry as a p(2  1) type structure rotated by 458 as shown in the top panel of Fig. 1. The thermal stability of the ordered phase structure can be partially monitored through an examination of the temperature dependence of the polar (u) and azimuthal (w) angles distributions. The polar angle (u) distribution is shown in Fig. 2, where nitrogen molecules tilted in average directions by the same polar angle of 538 with respect to the surface normal. As the temperature is

Fig. 2. Polar angle distribution is plotted for temperatures 1 K, 20 K and 25 K. At 1 K the distribution is symmetric and centered at u = 538.

Fig. 3. Azimuthal angle (w) distributions of the N2 molecules adsorbed on LiF surface, at temperatures 1 K, 20 K and 25 K, show the progress of the transition from an ordered to disordered phase.

2977

increased to around T = 23 K, no change in the position of the peak of the polar (tilt) angle distribution was observed. The w distribution (as shown in Fig. 3), indicates that at 1 K, there are two sharp peaks at w = 08 and 908. This indicates that the nitrogen molecules have the same azimuthal orientation along diagonals, with diagonals alternating their orientation. These peaks still exist up to 20 K even they have been thermally broadened. Beyond 20 K the azimuthal angle distribution is found to have four peaks at w = 908pffiffiffiandpwffiffiffi= 1808. Together these distributions indicate that the pð2 2  2Þ R 45 structure persists up to 20 K and it is clearly that the molecules start to lose their preferred azimuthal orientations at T > 20 K. These results show that the conversion of the ordered to the disordered phase is occurred at transition temperature that is tentatively assigned to the range from 20 K to 25 K, but our simulation results of the variation of the average molecular energy pffiffiffi of pffiffiffiN2 in adsorbed phase with temperature shows that the pð2 2   2Þ R 45 ordered phase converts to disordered phase starting at T = 23 K (see Fig. 4). Beyond 25 K, the system enters another new phase where an empty area in the layer is formed as shown in Fig. 5. This means that the layer is still willing to have more molecules at temperatures above 23 K where the ordered commensurate structure may be replaced by a disordered structure at higher density. This is unusual situation since denser structures are more ordered due to the

Fig. 4. The variation of the average molecular energy of the adsorbed nitrogen plotted for temperatures up to 35 K.

Fig. 5. The LiF(0 0 1) surface covered with 72 N2 molecules at 27 K. The nitrogen atoms are shown as gray circles. The (+) symbol represents a Li+ ion, whereas the () symbol represents an F ion.

2978

A.K. Sallabi et al. / Applied Surface Science 256 (2010) 2974–2978

motion restrictions caused by the close packing. This kind of behavior has been only seen in the adsorption of CO2 molecules on MgO(0 0 1) surface [1]. 4. Conclusion The interaction potential that adopted in this work yields a reasonable result for the adsorption of a single nitrogen molecule over a LiF surface in terms of binding energy, molecular orientation and agrees well with the available experimental findings. Therefore it was applied here to study the structural phase transition of N2/LiF(0 0 1) system at temperature range of 1–35 K using a Monte Carlo simulation method. The Monte Carlo simulations show that nitrogen molecules pffiffiffi form pffiffiffi an order structure on LiF(0 0 1) surface with a pð2 2  2Þ R 45 unit cell symmetry. This structure undergoes a phase transition into a disordered phase. It is found that above 23 K the molecules loss their azimuthal order but with residual short-range order. We observe no change in the position of the peak of the polar (tilt) angle distribution above the transition temperature and so this transition is purely of the order–disorder type. Acknowledgment DBJ wishes to thank the Natural Sciences and Engineering Research Council (Canada) for financial support in the form of a discovery grant.

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]

C.D. Daub, G.N. Patey, D.B. Jack, A.K. Sallabi, J. Chem. Phys. 124 (2006) 114706. A.K. Sallabi, D.B. Jack, Phys. Rev. B 61 (2000) R4841. A.K. Sallabi, D.B. Jack, J. Chem. Phys. 112 (2000) 5133. V.J. Barclay, D.B. Jack, J.C. Polanyi, Y. Zeiri, J. Chem. Phys. 97 (1992) 9458. J.C. Polanyi, R.J. Williams, S.F. O’Shea, J. Chem. Phys. 94 (1991) 978. J.V. Jose´, L.P. Kadanoff, S. Kirkpatrick, D.R. Nelson, Phys. Rev. B 16 (1977) 1217. G.E. Ewing, Int. Rev. Phys. Chem. 10 (1991) 391. J. Heidberg, E. Kampshoff, M. Suhren, J. Chem. Phys. 95 (1991) 9408. D. Schmicker, J.P. Toennies, R. Vollmer, H. Weiss, J. Chem. Phys. 95 (1991) 9412. P.N.M. Hoang, S. Picaud, C. Girardet, A.W. Meredith, J. Chem. Phys. 105 (1996) 8453. G. Dash, Films on Solid Surfaces, Academic, New York, 1975; G. Dash, Contemp. Phys. 30 (1989) 89. A. Lubezky, Y. Kozirovski, M. Folman, J. Phys. Chem. 97 (1993) 1050. K. Hermansson, M. Alfredsson, Surf. Sci. 411 (1998) 23. P.N.M. Hoang, S. Picaud, C. Girardet, J. Chem. Phys. 105 (1996) 8453. C. Graham, J. Pierrus, R.E. Raab, Mol. Phys. 67 (1989) 939. H. Hettema, P.E.S. Wormer, A.J. Thakkar, Mol. Phys. 80 (1993) 533. W.J. Meath, A. Kumar, Int. J. Quantum Chem. S24 (1990) 501. A. Kumar, W.J. Meath, Theor. Chim. Acta 54 (1992) 131. G.D. Zeiss, W.J. Meath, Mol. Phys. 33 (1977) 1155. C. Hattig, B.A. Hess, J. Chem. Phys. 105 (1996) 9948. T.L. Gilbert, J. Chem. Phys. 49 (1992) 2640; T.L. Gilbert, O.C. Simpson, M.A. Williamson, J. Chem. Phys. 63 (1975) 4061. F.T. Smith, Phys. Rev. A 5 (1972) 1708. W.R. Busing, WMIN: A Computer Program to Model Molecules and Crystals in terms of Potential Energy Functions, ORNL-5747 (available through QCPE). M.A. Sabri, MSc. Thesis, Concordia University, Montreal-Canada, 1996. K.T. Tang, J.P. Toennies, J. Chem. Phys. 80 (1984) 3726. K.T. Tang, J.P. Toennies, Z. Phys. D 1 (1986) 91. C. Douketis, G. Scoles, S. Marchetti, M. Zen, A. Thakkar, J. Chem. Phys. 76 (1982) 3057.