Structures of D 2 Layers on MgO(001) ∗

Monte Carlo (MC) simulations of D 2 molecules on the MgO(001) surface are reported and show that a series of interesting structures form with increasing coverage, viz. p (2 × 2) → p (4 × 2) → p (6 × 2), with coverages Θ = 0.5, 0.75, and 0.83 respectively, and are stable up to 13 K. The p (2 × 2) structures contain two D 2 molecules per unit cell, with each molecule lying parallel to the plane of the surface ( θ = 90 ◦ ) directly above every other Mg 2+ site. The molecules adopt a “T” conﬁguration with respect to their nearest neighbors. The p (4 × 2) and p (6 × 2) structures, have two kinds of adsorption sites: a parallel site, as in the case of p (2 × 2), and a tilted site, where the D 2 molecules sit between cationic and anionic sites with the molecular axis directed towards the anionic site, with θ ≈ 60 ◦ . These structures are consistent with recent neutron scattering results in terms of coverage and stability, but disagree in terms of symmetry; the neutron scattering work found “ c ” type structures whereas the MC simulations (without quantum considerations) yield a “ p ” type structures. To reconcile the results of the simulations and experiments, the quantum mechanical rotational motion of the adsorbed D 2 molecules was studied using perturbation theory. These calculations show that the adsorbed D 2 molecules are azimuthally delocalized and hence the structures are indeed “ c ” type rather than “ p ” type. [DOI: 10.1380/ejssnt.2009.207]


I. INTRODUCTION
The adsorption of D 2 molecules on the MgO(001) surface has been of interest for some time and has come under investigation using neutron scattering technique [1]. Although one of the simplest adsorbate-substrate systems imaginable for the study of gas-solid interactions, it has the additional possibility of exhibiting quantum effects, owing to the light mass of the adsorbate. Early studies using neutron scattering [1], adsorption-isotherm and heat capacity measurements [2], and nuclear magnetic resonance [3] have examined the structure of submonolayer films of H 2 , HD, and D 2 on MgO(001) surface. The neutron scattering study found evidence of a sequence of phase transitions c(2 × 2) → c(4 × 2) → c (6 × 2) in the D 2 /MgO system. This sequence was recently observed for all three adsorbate systems (H 2 , HD, and D 2 ) and information about the adsorbate vibration frequencies extracted using helium atom scattering (HAS) [4]. The details of the adsorbate structure are not directly accessible by experimental means. So theoretical and computational methods must be used to clarify these information since the details of the molecular positions are needed to complete our understanding of these structures and will also set the stage for a discussion of quantum effects, namely, the role of molecular rotational states (ortho and para) in the adsorbed phase.
Previous Monte Carlo (MC) work on molecular adsorption on surfaces has yielded good agreement with experimental findings and hence the same methods were used to determine the details of these structures. Here we present the results of Monte Carlo (MC) simulations combined with perturbation theory (PT) calculations of molecular D 2 physisorbed on a MgO(001) surface.

II. INTERACTION POTENTIAL
To study the structure of a D 2 monolayer on MgO(001) theoretically and computationally, it is necessary to construct a reasonable potentials of molecule-molecule and molecule-surface interactions. Note here the interaction potential between deuterium molecules in its ground state is approximately determined, whereas there are no potentials in any other rotational levels (J > 0). Thus classical pair potentials provide a reasonable route to simulate a large numbers of molecules over the surface. The methods and potentials that applied in this work was explained in details in the simulations of the CO/MgO(001) and H 2 /LiF(001) systems [5,6]. Therefore, only an outline of the interaction potential calculations will be presented below.

A. Molecule-molecule (D2-D2) interactions
Since molecular D 2 has only a quadrupole moment, the D 2 molecule is represented as atomic point dipoles of strength 0.434 D placed at each atomic site. These atomic dipoles are directed along the D 2 bond axis away from the center of mass to reproduce the experimental molecular quadrupole moment of +0.6423 D·Å and yield a zero molecular dipole moment [7]. Therefore the electrostatic interaction among these molecules over the surface was calculated as two body summation of inter-atomic dipoledipole interaction, whereas, the inter-molecule dispersionrepulsion interactions was calculated using Buckingham formula, U (r) = A exp(−β · r) − C 6 /r 6 , where r is the inter-nuclear distance between deuterium atoms, C 6 is the London dispersion coefficient, and the parameters A and β are the Born-Mayer parameters. In this work, the repulsion-dispersion interaction between D 2 molecules is assumed to be the same as that for H 2 -H 2 interaction as assumed elsewhere [8]. Consequently, we are using the repulsion-dispersion parameters of molecular H 2 that estimated from data obtained from ab initio studies of H 2 dimers [9,10] and are found to be A D = 980.89 kcal/mol, β D = 3.212Å −1 , and C 6 = 6.51Å 6 ·kcal/mol.

B. Molecule-surface (D2-MgO) interactions
The total interaction between a D 2 molecule and the substrate is assumed to be a pair-wise sum of contributions from electrostatic, dispersion and repulsion interactions. The electrostatic contribution represents the interaction between the D 2 atomic dipoles and the surface electric field plus the induction energy. The repulsion-dispersion interactions between the atoms of a D 2 molecule and the ions of the substrate are modeled using Tang-Toennies potential [11]. The repulsiondispersion parameters of a D atom interacting with Mg 2+ /O 2− ions were developed using the same methods as applied in Ref. [6] and are presented in Table I. Our model considers D 2 as a rigid rotor of bond length 0.74Åinteracting with a stationary (001) MgO lattice, whose the distance between nearest like ions in the surface is 2.97Å. A two-dimensional Fourier series was used to represent the molecule-surface interaction potentials as previously applied successfully to the HBr/LiF(001) system [12,13]. These assumptions greatly simplify the Monte Carlo energy calculations and the molecule is now a rigid body sitting in a periodic potential. Over the surface, the orientation of the deuterium molecule was defined by a polar angle θ (tilt with respect to the surface normal) and an azimuthal angle ϕ (angle in the plane of the surface).

III. SIMULATION METHOD
Our computational technique employs the familiar canonical Metropolis Monte Carlo (MC) method [14], where randomly selected molecules are allowed to change their positions and orientations in accord with the Boltzmann criteria. In our simulations, each cycle consisted of every molecule attempting a move. In this way, many molecular classical configurations are sampled and the potential energy of each one is calculated, the various configurations are then weighted according to the Boltzmann distribution and obey classical statistics.

IV. SIMULATION RESULTS
The atom-ion parameters that are presented in Table I has been tested by studying the adsorption of a single D 2 molecule on MgO(001) using the Monte Carlo simulations. On average, the molecule prefers to sit flat (θ = 90 • ) on the top of the Mg 2+ site at a height of 2.76Å where the molecular axis of D 2 is oriented towards the closest neighboring O 2− sites (ϕ = ±45 • ), with a potential well depth of −0.605 kcal/mol. The perpendicular motion of a D 2 molecule to the surface is expected to be quantized with an appreciable zero point energy due to its light mass. Therefore, the surface potential, described above, was fit to a 1-D harmonic oscillator model, and the zero point energy was estimated to be 0.065 kcal/mol (z-motion frequency) which yields a binding energy of −0.54 kcal/mol. Note here there is no available experimental data for the binding energy of D 2 /MgO system, but we can compare it with that of H 2 molecule since the intermolecular interaction between those species are approximately identical [8]. Therefore the calculated value of the binding energy of a single D 2 molecule over MgO surface falls within the range estimated by Skofronick (−0.461 to −0.577 kcal/mol) [4] and still in good agreement with other theoretical calculations of −0.606 kcal/mol [15]. Our finding that D 2 , at low coverage, prefers to sit at the top of Mg 2+ site is consistent with proposals based on the neutron scattering results [1] and hence our potential model seems to be a reasonable one to study the structures of D 2 layers on a MgO surface.
The neutron scattering results [1] show that at a coverage of 0.5, an ordered commensurate structure with c(2 × 2) symmetry appears. Degenhardt's group proposed a unit cell structure with deuterium molecules over every other Mg 2+ sites arrayed in a checkerboard fashion [1]. This structure and others were examined using the Metropolis Monte Carlo method, in which a patch of surface with 144 Mg 2+ sites and 72 D 2 molecules are used to yield the correct coverage. These simulations were run for 60 kcycles with the first 10 kcycles thrown away. Regardless of the initial configurations, a unique final structure was obtained where the molecules are sitting flat above Mg 2+ sites and arranged in a unit cell of a p(2×2) symmetry as illustrated in Fig. 1. This structure remained stable up to 13 K. Beyond that temperature, molecules migrated from their adsorption sites without diffusion along the surface and some of them started completely to desorb from the surface. Thus the p(2×2) structure is unstable beyond 13 K.
At densities above 50% monolayers, neutron scattering technique found two other phases. The c(4 × 2) phase has an estimated coverage of 0.75, where the molecules are equally spaced along the cationic channels of the sur-  face [1]. This proposed configuration and many other configurations with different molecular orientations were examined using the MC simulation. In spite of the chosen initial configurations, the simulations all evolved to a single final configuration with a p(4 × 2) symmetry as shown in Fig. 1. The unit cell consists of two kinds of adsorption sites; a parallel site, as in the case of p(2 × 2), and a tilted site, where the center of mass of D 2 molecules lie between the cationic and anionic sites with their molecular axes oriented towards the nearest O 2− sites with a polar angle of 57 • . At T = 13 K, the molecules start to leave their adsorption sites causing disorder within the layer and hence the p(4 × 2) layer are stable up to 12 K. At a coverage of 0.83, the neutron scattering results show the existence of a third monolayer where the molecules arranged in a unit cell of c(6 × 2) symmetry and consists of two types of adsorption sites [1]. This proposed structure and other configurations with different molecular orientations were examined using the MC simulations. Our results confirmed the existence of a third monolayer with a unit cell of a p(6×2) symmetry as shown in Fig. 1. This monolayer contains two kinds of adsorption sites, parallel and tilted, similar to that found in case of p(4×2). Within a unit cell there are eight tilted molecules and two flat molecules for a total of ten molecules. This monolayer was found to be stable up to 12 K, where the molecules remain in their adsorption sites and exhibit in average the same molecular orientations over the surface. At T > 12 K, the tilted molecules migrate from their lattice sites and would desorb from the surface plane. Thus the p(6 × 2) layer persists up to 12 K.
The results obtained by using the Metropolis Monte Carlo method match the experimental results in terms of coverage and stability but disagree in terms of symmetry. In particular, p-type structures are observed rather than c-type due to the localization of the molecular axes along either the ϕ = 45 • or ϕ = 135 • azimuths. However, c-type structure can be obtained if the molecular axis of deuterium is azimuthally delocalized or freely rotating. This cannot be observed in our simulations since the MC method uses classical statistics and is unable to yield stable c-type structures. Therefore the rotational motion of a molecular D 2 in adsorbed phase was examined using perturbation theory to see if delocalization would occur and how adsorption site preference depends on the rotational state of D 2 (ortho-vs. para-species).

V. ROTATIONAL STATE CALCULATIONS
For a D 2 molecule sitting at fixed position directly at the top of Mg 2+ or O 2− sites and whose molecular axis can be defined in terms of θ andϕ, the molecule-surface rotational potential has the following formula, The above formula reflects the fourfold symmetry of the surface as well as the rotational symmetry of the homonuclear symmetry of deuterium molecule. The A n coefficients were determined using the square fit program from Maple V. 7 and the values are listed in Table II.
As done in Ref. [6], the wave functions of D 2 in different rotational state were estimated using time-independent perturbation theory to the first order correction level and may be calculated using the following formula; Over the Mg 2+ site, the angular distribution of the ground state ortho-D 2 and helicoptering para-D 2 (J = 1, m J = ±1) states have similar behavior with a maximum value when the molecule lies flat over the surface (θ = 90 • ) and a minimum when it lies vertical to the surface plane (θ = 0 • ) as shown in Fig. 2. Consequently, ortho-D 2 is highly squashed over the Mg 2+ site, in sharp contrast to its uniform angular distribution in gas phase. At θ = 90 • , the azimuthal distribution is almost uniform with a small ripple that reflects the fourfold symmetry of the surface. In other words, the orientation of the molecular axis is azimuthally delocalized with a residual preference in the ϕ = ±45 • , ±135 • directions. For cartwheeling para-D 2 (J = 1, m J = 0), the angular distributions has a maximum when the molecular axis adopts the surface normal orientation, and a minimum when its molecular axis makes a right angle with the surface normal (results not shown). Over the O 2− site, the angular distribution of the ground ortho-D 2 and cartwheeling para-D 2 (J = 1, m J = 0) states are similar. The angular distribution has a minimum value when the molecule locates flat over the surface and a maximum value when its axis lies parallel to the surface normal (see Fig. 3). Compared to the results obtained above the Mg 2+ site, the angular distribution of the ground ortho-D 2 is completely reversed; ortho-D 2 is elongated over the F − site. With respect to the az- imuthal angle, a uniform angular distribution is obtained with small ripple that reflects the fourfold symmetry of the surface. For helicoptering para-D 2 (J = 1, m J = ±1), the angular distributions also resemble their gas phase counterparts (results not shown).
Regardless of a molecule's rotational state or adsorption site, the angular distributions show that the molecular axis is azimuthally delocalized. Applying this to the results of the MC simulations, we note that the molecules located at the center and corners of the unit cells now lose their preferred azimuthal orientation and become equivalent. The p-type structures thus become c-type structures, and the MC results are reconciled the experimental findings.
Perturbation theory also yields the rotational energy levels, namely where the spherical harmonics have been replaced by the usual separable functions as This formula is used to calculate, for any given values of the quantum number J and m J , the rotational energy levels above either a cation or anion. Note here that at low temperature (T < 20 K), the J = 0 and J = 1 rotational levels are, practically speaking, the only ones populated according to the equilibrium Boltzmann distribution and hence we are only interested in calculating the rotational energy for these two levels. The energies of molecular D 2 in different rotational levels above Mg 2+ site are calculated using above formula and are listed in Table III. In general, the adsorption onto the surface removes the degeneracy present for a given rotational state J. Each energy level is split into an ordered sequence of levels. Within a set of splittings, the |m| = J states always have an energy lower than the corresponding gas phase value as well as the lowest energy in the set, while the m = 0 state always has the highest energy and when J = 0, these states have an energy higher than the gas phase. These observations indicate that the |m| = J states are allowed to adsorb on the surface (as well as the ortho J = 0 state), whereas the m = 0 states cannot adsorb.
The energy states of a D 2 molecule adsorbed above an O 2− site are also listed in Table III. Again it is found that each rotational energy level J splits into an ordered sequence of levels. For a given set of splittings, the adsorbed m = 0 states have a lower energy than the corresponding gas phase value, while the |m| = J states have the highest energy in the set, and always have energies higher than the gas phase. Consequently, the m = 0 states are allowed to adsorb over the anionic site whereas the |m| = J states  [4,6], even the mass of D 2 is increased by a factor of 2 with respect to H 2 , our PT calculations show that ortho-D 2 and para-H 2 species exhibit the same behaviour in adsorbed phase, in terms of adsorption site, cylindrical distribution and adsorbate structures. In particular, the helicoptering states of para-D 2 and ortho-H 2 species prefer to sit flat at the top of the cationic site whereas the cartwheeling states of para-D 2 and ortho-H 2 species prefer to sit along the surface normal at the top of the anionic site. On MgO surface, the adsorbate structures of H 2 and its isotope D 2 should be the same since the three kinds of commensurate mono-layers of D 2 and H 2 , c(2 × 2) → c(4 × 2) → c(6 × 2) with coverages Θ = 0.5, 0.75, and 0.83 respectively, are observed experimentally using different techniques [1,4].

VI. CONCLUSION
Through the combination of Monte Carlo simulations and perturbation theory calculations of rotational states we have shown that the sequence of structures, c(n × 2) previously observed experimentally is indeed possible. Furthermore, the details of these structures have been determined. The c(2 × 2) structure consists of an array of molecules covering every other Mg 2+ site of the surface in a checkerboard pattern, with quantum mechanical delocalisation of the molecular axes eliminating azimuthal differences between molecules. Specifically, the ortho and helicoptering para states are allowed to adsorb here. The c(4×2) structure consists of two kinds of adsorption sites. One third of the molecules adsorb directly over Mg 2+ ions with a preference for a horizontal orientation for the molecular axes (ortho or helicoptering para-states), while the remaining two thirds adsorb near, but are offset from, the O 2− ions with orientations that prefer a tilt from the surface normal. In terms of rotational states these are thought to be a mix of cartwheeling and helicoptering para-states or possibly skewed ortho-states. The tilted molecules sit 0.5Åfurther from the surface than the horizontal molecules. The c(6 × 2) structure is an extension of the c(4 × 2) structure with only 1/5 of the molecules adopting a horizontal orientation; the rest are tilted near O 2− ions.