Oxidative etching mechanism of the diamond (100) surface

We performed density functional theory calculations with van der Waals corrections to elucidate diamond oxidation mechanism on the atomic-level which could lead to insights that will advance the improvement of nascent nanofabrication technologies. We developed a comprehensive theory of oxidative etching of the diamond (100) surface, from the adsorption of gas phase O2, including details of metastable adsorption states, intersystem crossing, and induced surface dereconstruction, to the desorption of CO and CO2, complete etching of the top surface layer and its subsequent stabilization. Oxygen adsorption induce C-dimer bond breaking through the formation of carbonyl structures at low surface coverage. At high surface coverage, adjacent carbonyls can transform into the more stable ether chain with small energy barrier requirement. A mix of carbonyl and ether will form in surfaces with defects, and several models of possible structures have been presented. CO desorption creates a point defect that serves as nucleation site for the preferred etching direction along [011], perpendicular to the top layer C-dimer bonds. We obtained a wide range of desorption activation energies, which depend on the existence of point defects and reconstruction of surfaces. C-dimer formation and oxygen adsorption stabilize vacancies and single-atomic-layer-deep trough. © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Diamond possesses a unique combination of superlative properties that makes it an attractive material for several novel applications. Of all semiconductors, diamond has the highest dielectric breakdown field, saturated electron drift velocity, and thermal conductivity which make it an ideal material for next-generation electronics that will outperform current systems in terms of operating frequency and power handling capacity [1]. In addition, diamond's wide band gap, optical transparency, chemical inertness, hardness, and low thermal expansion have led to a number of specialized photonics applications including spacecraft window, optoelectronic microlens, photodetector, high energy particle detector, UV and infrared LED, and Raman laser [2]. Moreover, the low density and coefficient of friction of single crystal diamond makes it suitable for micromechanical and microelectromechanical applications [3,4]. Diamond has color centers that functions as bright single photon source (SPS) at room temperatures and can be potentially implemented in quantum metrology and quantum information science [5e8]. Furthermore, diamond is non-toxic and biocompatible. Currently, several applications of diamond in biotechnology is being explored, such as drug delivery, gene transfection, fluorescence imaging, and magnetic resonance imaging [9,10].
Extreme hardness and chemical stability are the key properties that give diamond distinctive advantage over other materials. However, these properties make fabrication of diamond nanostructures exceptionally challenging. To utilize diamond's full potential, a complete fabrication toolkit similar to the ones available for silicon is required [2]. A number of diamond nanofabrication techniques based on oxidative etching are being developed. The most widely used is oxygen plasma reactive ion etching (RIE) process. Several diamond devices were fabricated using RIE including vertical-type metal-oxide semiconductor fieldeffect transistors (MOSFET) [11], vertically aligned 1-dimensional diamond nanostructures [12e14], optical nanocavities [15], nanocrystalline resonators [16], and micro-mechanical components [3]. Fine-tuning of the ratio of O 2 to other gases like Ar in inductively coupled plasma RIE processes leads to etched surfaces with higher quality [17]. Nevertheless, plasma etching has some drawbacks including low diamond-mask selectivity and the potential to damage diamond devices that will cause performance deterioration [18]. An alternative method is thermochemical etching, which implements the observed catalytic function of transition metals such as Fe and Ni on diamond tool wear [19,20]. Morofushi et al. thermochemically etched single crystal diamond with Ni mask under high-temperature air and reported that diamond etching was caused either by carbon diffusion into nickel or oxidationreduction reaction between nickel oxide and carbon. Significantly lower etching rates were recorded when diamond is etched under nitrogen gas highlighting the role of oxygen in accelerating the etching rates [21]. Nagai et al. performed similar thermochemical etching using high-temperature water vapor and reported anisotropy in etching of low-index diamond surfaces. The (100) surface under Ni film was etched with inverted pyramid and truncated pyramid-shaped trenches with {111} side walls functioning as stopping planes. On the other hand, the (110) surface was etched with U-shaped diamond trenches and vertical {111} sidewalls. In contrast to the two planes, the (111) surface was atomically flattened, albeit no trench was formed [18,22].
Oxygen plays a significant role in both plasma and thermochemical etching, thus detailed understanding of the interaction of oxygen with diamond surface will lead to the improvement of the quality of these fabrication techniques. The (100) surface of diamond, C(100), is one of the dominant surfaces produced during chemical vapor deposition (CVD) synthesis and is the most technologically relevant surface, specifically for microelectromechanical systems (MEMS) and power electronics device fabrication [11,23e28]. The adsorption of oxygen on the C(100) surface has been studied extensively both experimentally and theoretically. Electron diffraction measurements have revealed that oxygen adsorption induced breaking of surface carbon dimer bond, leading to the dereconstruction of the C(100) surface [29,30]. Electron energy loss spectra (EELS) and high-resolution x-ray photoelectron spectroscopy (XPS) of the C(100) surface exposed to atomic and molecular oxygen have demonstrated that the C(100) surface is terminated by oxygen atoms on top and bridge sites, forming carbonyl and ether groups, respectively, with coverage of up to a monolayer [31,32]. Theoretical studies have shown that while both top site and bridge site are structurally stable at monolayer coverage, the latter site is energetically preferred [23,33e37]. Recently, Lu et al. proposed a methoxyacetone surface structure to account for the observed co-existence of carbonyl and ether [37].
In addition, several studies were performed to investigate diamond etching using dry oxygen at high temperatures. C(100) surfaces etched with oxygen were generally flat with very low steps whose sides are along the [110] direction. In contrast, the etched C(111) surfaces were morphologically rough [38]. Moreover, atomic force microscopy (AFM) images of C(100) have shown that the surface roughness has no appreciable change after removal of 480 nm of carbon [32]. Thermal oxidation experiment has also shown that the C(100) surface is more oxidation resistant than the C(111) surface [39]. This has been confirmed through a molecular dynamics simulation of energetic atomic oxygen collision with lowindex diamond surfaces [40]. Analysis of thermal desorption spectra (TDS) revealed that the oxidative etching of diamond surface is mainly by desorption of CO although much smaller amounts of desorbed CO 2 were also detected [31,41]. The most puzzling feature of the TDS data is that the peak of CO desorption was extremely wide. From the peak temperature of TDS, Frenklach et al. estimated the CO desorption energy to be 45.0 kcal/mol (1.95 eV). In addition, they performed semiempirical quantum chemical potential energy calculations to investigate the CO desorption mechanism. The largest potential energy barrier is 38.4 kcal/mol (1.67 eV) which corresponds to the breaking of the final CeCO bond. While the calculated potential energy barrier agreed reasonably with their estimate from TDS data, their simulation was not able to explain the observed wide desorption peak [33].
At present, a theoretical study that investigates the entire oxidation process, starting from the adsorption of gas phase O 2 up to the desorption of CO and etching of the surface, has never been done. This kind of study is necessary to develop a comprehensive knowledge of diamond oxidation mechanism. For this reason, van der Waals-corrected density functional theory simulations of the oxidation of the C(100) surface have been performed. The intersystem crossing, O 2 dissociation, surface dereconstruction, differential adsorption energies, desorption of CO and CO 2 , etching of the entire top surface layer, and etched surface stabilization were all studied rigorously. The results of our calculations were analyzed and compared with available experimental data to gain detailed understanding of the interaction of the C(100) surface with oxygen.

Model of the diamond (100) surface
The calculated lattice constant for bulk diamond is 3:572 A which agrees very well with experimental value of 3:567095±1:7 Â 10 À5 A [50]. The C(100) surface was modeled as a slab of 8 carbon layers with a vacuum region of 35 a 0 (18:5 A) which is about 13 carbon-layers thick. The top layer C atoms were dimerized to conform to the 2 Â 1 reconstruction observed in electron diffraction measurements [31,51]. The bottom layer was reconstructed in a similar manner and terminated with hydrogen atoms. Earlier theoretical works on the electron affinity and surface conductivity of C(100) with various surface functionalizations have shown that thicker slab is necessary to describe the bulk electronic structure [52,53]. For this reason, we used a thicker slab with 20 layers of carbon atoms for band structure calculations. Comparison of the band structures for 8, 12, 16, and 20-layer models is shown in the Supplementary Document (Fig. S1-S2). Calculations have also shown that the adsorption energies require smaller number of layers to converge (Table S1-S4). Optimizing the geometry of a nonreconstructed C(100) surface will lead to the same 2Â 1 reconstruction, suggesting that the process happens spontaneously. The calculated C-dimer bond length is 1:383 A. This agrees well with earlier theoretical studies that reported dimer lengths between 1:36 À 1:38 A [54e57]. For all adsorption and oxidation simulations, the top five carbon layers were allowed to relax, while the bottom three carbon layers and H terminating layer were fixed.
A 2 Â 3 supercell with 12 Â 8 k-point mesh was used in oxygen adsorption simulations to limit the interaction of oxygen species with adjacent repeat unit. For differential adsorption energy calculations, a larger 2 Â 4 supercell with 12 Â 6 k-point mesh was used to model various adsorption coverages. Similarly, this larger supercell was used for oxidation simulations, as it was found to be necessary for highly detailed modelling of the CO desorption mechanism (Fig. 1).

Adsorption energies
The adsorption energy E ads is defined with respect to the energies of clean 2 Â 1 reconstructed surface E surf and isolated gasphase triplet O 2 E O2 : where E sys is the energy of the system of diamond surface and adsorbed O atoms. In this definition, the positive and negative values of E ads mean repulsive and attractive interaction, respectively. The differential adsorption energy of O atom on the C(100) surface with respect to gas-phase O 2 molecule DE O = O 2 is given by: where E surf þnO2 and E surf þðnÀ1ÞO2 are the total energies of system with n and n À 1 number of O 2 adsorbed, and E O2 is the total energy of gas phase O 2 . With respect to gas-phase O atom with total energy E O , the differential adsorption energy DE O =O is given by:

Crystal orbital overlap population analysis
The bonding mechanism of O and CO on diamond surfaces were analyzed using the population analysis scheme applicable to DFT calculations with plane wave basis sets [58,59]. In this scheme, the total system is divided in two subsystems A and B corresponding to CO and diamond surface respectively. A single-electron wavefunction j i is expressed in terms of the sum of the minimal basis sets N A and N B for systems A and B where S j;k is the overlap integral between wave functions c A j and c B k .

Vibrational properties
The vibrational properties were calculated using finitedifference harmonic approximation by diagonalizing the Hessian matrices to obtain vibrational eigenvalues and eigenvectors. The Hessian matrices are numerically estimated by initially displacing the O atoms and top four C layers by 0:05 a 0 . The displacements of each atom were updated until the energy eigenvalues converge. Two methods were implemented to produce vibrational spectra from the calculated eigenvalues. First, Lorentzian broadening with full-width-half-maximum (FWHM) of 15 meV were used to match the resolution of reported EELS data. For the second method, we calculated the intensity of vibration EELS peaks using dipole scattering mechanism [60]. In this scheme, the relative intensities of the vibrational modes are assumed to be proportional to the square of the dynamic dipole moment dm=dQ I loss I elastic f dm dQ 2 (8) The dynamic dipole moment was calculated from the difference between the work function on the two sides of the simulation slab [61e63]. Effective screening medium (ESM) based on Green's function method was implemented to eliminate the non-physical dipole and multipole interactions with image slabs [64,65]. The Bohr radius, vacuum permittivity, and normal frequency are given by a 0 , ε 0 , and u s , respectively. This calculation method involves the experimental parameters of HREEL spectroscopy, namely the primary energy E I , incident angle of electron beam q I , and acceptance angle of spectrometer q C , which defines b q C ¼ q C =q E , where q E ¼ Zu s =2E I . The experimental parameters affect the calculated intensity of the spectra but do not alter the main conclusion of this work. The coverage of adsorbates is given by n s , while the function Finally, the intensity data were broadened using Gaussian with FWHM of 8 meV.

Results and discussions
3.1. Adsorption of gas-phase O 2 on the 2 Â 1 reconstructed diamond (100) surface In this section, we simulated the adsorption of gas-phase O 2 on top and bridge sites of the clean, reconstructed diamond (100) surface, C(100)-(2 Â1). We start by performing test calculations to guide our choice of method that will account for dispersion relation. Fig. 2 shows the potential energy curve of singlet O 2 adsorption on the top site of the C(100)-(2 Â1) surface calculated with constrained coordinates using the PBE exchange correlation functional, and GGA functionals with semiempirical dispersion correction (PBE þ D2) [46], and nonlocal correlation corrections (optb86b-vdW and rev-vdW-DF2) [66,67] to account for dispersion forces. The dispersion corrections enhanced the adsorption energies by À0:241 eV, À0:354 eV, and À0:284 eV for PBE þ D2, optb86b-vdW, and rev-vdW-DF2, respectively. On the other hand, the adsorption distance, d predicted by optb86b-vdW is identical to PBE, while the d for PBE þ D2 and rev-vdW-DF2 slightly decreased and increased, respectively by 0:001 A. The test calculations have shown that dispersion correction is significant for our system. We decided to use the PBE þ D2 functional as it gives the best compromise for accuracy and computational cost.
Next, we begin with the simulation of the initial oxidation process. Gas phase O 2 is a spin-polarized triplet ( 3 S À g ) molecule in its ground state. As O 2 approaches surfaces such as silicon and carbon nanotubes, intersystem crossing (ISC) to the more reactive singlet state ( 1 D g ) occurs [68e70]. Therefore, spin polarization is necessary to properly describe the initial oxidation process of diamond. In our calculations, triplet state was obtained by fixing the difference between the number of spin-up and spin-down electrons to two. As for the singlet state, spin restricted calculation was performed to obtain a closed shell configuration (spin-restricted) which is necessary to prevent spin contamination [69,70]. Our calculated triplet-singlet excitation energy is 1:17 eV which agrees reasonably well with the experimental value of 0:98 eV [71].
The calculated adsorption reaction paths of triplet and singlet O 2 on the top site are shown in Fig. 3 using red and blue curves, respectively. A metastable physisorbed state is observed in the reaction path of triplet O 2 with adsorption energy of À0:78 eV. From the metastable state, the system must overcome a 0:55 À eV energy barrier to reach the final molecular adsorption state (E ads ¼ À1:67 eV). In comparison, the reaction path of singlet O 2 is barrier-less with more stable final state (E ads ¼ À 2:67 eV).
Triplet -singlet O 2 conversion is spin-forbidden but it can be induced by spin orbit coupling within the O 2 /diamond interface. The probability of this transition is low because spin-orbit interaction is weak in light elements [71]. Triplet O 2 is first trapped in a chemisorption well, as in the case of initial oxidation of Si (100) surface. Although the probability of spin transition is low, as the adsobed O 2 undergoes molecular vibration, it passes to the ISC point many times until it switches to the singlet reaction path [70]. To gain further insight, we simulated the ISC process. The transition from triplet state to singlet state occurs at the ISC point, where their energies are degenerate. The geometry at the ISC point was initially estimated at the crossing of the triplet O 2 and singlet O 2 reaction paths. Subsequently, the reaction path was calculated with replicas starting from the initial state (IS) up to the ISC treated with tripletconstrained multiplicity, while the remaining replicas up to the final state (FS) were simulated using closed shell model. The reaction path of O 2 that undergone spin transition is shown using black curve. Spin transition paves an easier path to the final adsorption state as the energy barrier is lowered to 0:27 eV. Nonetheless, this value suggests that spin transition is still an activated process. Thus, at initial oxidation, the oxygen molecules will be trapped in a metastable state. A part of the kinetic energy gained at finite temperature will be converted to rocking vibration of O 2 with one of its ends fixed at the surface. Even with low spin conversion probability, vibration in this manner will cause the  Similar calculations were performed for bridge site adsorption and the reaction paths are shown in Fig. 4. For the reaction path of triplet O 2 , after O 2 dissociation, a metastable chemisorbed state exists (E ads ¼ À 3:12 eV), where one O atom forms epoxide (cyclic ether) group with a C-dimer while the other bonds with the C atom of an adjacent C-dimer. The metastable state is separated from the final dissociative adsorption state (E ads ¼ À 3:35 eV) by a 0:78 À eV energy barrier. The reaction path of singlet O 2 is barrier-less and leads to a more stable configuration (E ads ¼ À 4:28 eV), similar to the reaction path on the top-site adsorption. ISC lowers the energy barrier of the dissociative adsorption to 0:35 eV.
Since the reaction path of triplet O 2 up to the metastable state is barrier-less for both top and bridge site adsorption, it is possible that both of these processes take place at initial oxidation. In addition, in both paths, the ISC takes place and the difference in the activation barrier height between the top site (0:27 eV) and the bridge site (0:35 eV) adsorption is small (0:08 eV). Therefore, it is likely that O 2 can be both adsorbed molecularly on the top site and dissociatively on the bridge site prior to diamond surface dereconstruction.

Oxygen-induced surface dereconstruction and oxygen adsorption at monolayer-coverage
Next, we investigated the oxygen-induced dereconstruction of C(100)-(2 Â1) surface as reported on electron diffraction experiments. This process involves breaking of C-dimer bonds and translating C atoms into the original bulk terminated positions. Fig. 5 shows the adsorption reaction path of O 2 from gas-phase (A) to molecularly adsorbed (D) state, through metastable (B) and spin transition (C) states, described in the previous section. Dereconstruction of the C-dimer with adsorbed O stabilizes the system by~2:4 eV (E). The calculated activation energy of the reaction is 0:81 eV which is much lesser than the strength of typical CeC bonds (3:6 eV) [72]. This proves that O 2 molecularly adsorbed in top position facilitates surface dereconstruction. As O 2 bond is broken, CeO bond order is increased and carbonyl group is formed. The energy of the system can be further lowered by translating the O atoms from the top site to the bridge site to form a chain of ether along the [011] (parallel to the dimer) direction (F). This reaction is expected to occur at higher coverages, where a row of carbonyl forms along the [011] direction. The formation of this ether chain requires a small activation energy of 0:13 eV. If carbonyls are not located at nearest neighbor to each other along this row, the reaction will not occur.
The adsorption of O 2 on the bridge site is shown in Fig. 6 (A-D). Fig. 6 (D-E) shows the C-dimer bond breaking reaction which transforms the epoxide group to ether group. This reaction is endothermic as it creates dangling bonds on the surface carbon atoms. These dangling bonds can be saturated with the adsorption of second oxygen molecule leading to the more energetically favorable chain of ether groups (D-F-G). The ISC takes place at the transition state (F). This chain structure is expected to form at higher surface coverages.
The differential adsorption energy DE O =O 2 of O atom on the C(100) surface with respect to gas-phase O 2 was calculated to estimate the likelihood of further oxygen molecule adsorption on oxygenated surfaces with various coverages (Fig. 8a). For coverages lesser than 1 ML, we assumed that O 2 dissociates and adsorbs on the bridge site and form ether chain. Therefore, the results reflect an upper-limit estimate of adsorption energies. For coverages beyond 1 ML, stable structures were obtained where O adsorbs on the backbond, bridging the firstand second-layer C atoms (Fig. 8c). For coverages up to 1 ML, the values range from À2:81 eV to À2:38 eV which predicts favorable adsorption. The small variation between these values suggests that the steric repulsion between adsorbed O on adjacent rows along [011] direction appear to be negligible. In addition, the surface dereconstruction on a particular row has minimal effect on the stability of adjacent rows. For coverages beyond 1 ML, DE O =O 2 is positive suggesting the difficulty of O adsorption on the backbond. This agrees with an experiment which has shown that the maximum O coverage for oxidized C(100) using dry O 2 is only up to a monolayer [32]. An earlier theoretical work has also shown that ether formation through backbond adsorption on second or third layer leads to large distortion of the crystal structure and predicts endothermic O 2   adsorption reaction [54]. However, in the case of C(100) oxidation using O atom, the possibility of reaching coverages beyond 1 ML could not be precluded. Fig. 8b  3.3. Electronic analysis of the diamond (100) surface with monolayer oxygen coverage Electronic analysis has been performed to gain insights on the bonding mechanism of O on the C(100) surface. As described in Section 2.2, we used a thicker 20-layer carbon slab model for band structure calculations to properly describe the bulk electronic states. The band structures of clean and oxygenated C(100) surfaces are shown in Fig. 9. The C(100)-(2 Â 1) surface is an indirect band gap semiconductor with isolated valence band maximum (VBM) and conduction band minimum (CBM) which correspond to p bonding highest occupied molecular orbital (HOMO) and p antibonding lowest unoccupied molecular orbital (LUMO) of surface carbon dimer, respectively (Fig. 9a). Adsorption of O on top site significantly narrows the energy gap to almost zero (Fig. 9b). The VBM corresponds to O lone pair orbital that has opposite phase with the s bonding orbitals between the first and second-layer C atoms, indicating the p antibonding character along OeC bond (Fig. 10a). In contrast, adsorption on the bridge site widens the energy gap and changes diamond (100) into a direct band-gap semiconductor (Fig. 9c). The VBM also corresponds to O lone pair and s orbitals but the antibonding interaction is expected to be weaker compared to that of top site adsorption because the distance between the O lone pair and s orbitals are farther (Fig. 10b).
Therefore, the VBM of the bridge site adsorption is more stable than the VBM of the top site adsorption.
To gain further insights on the bonding mechanism, we calculated the gross populations (GPOP) of O atomic orbitals and the corresponding overlap populations with respect to the bond between O and the surface, also known as the crystal orbital overlap population (COOP) (Fig. 11). The GPOP of 2p orbitals of O on the bridge site is more delocalized compared to the O on the top site, which suggests stronger hybridization in the bridge configuration. This can be confirmed with the COOP graph which shows wider bandwidth of bonding orbital overlap in the bridge site configuration compared to the top site configuration. The integral of COOP up to the Fermi level gives the total orbital overlap which scales directly with the strength of the bond. The total orbital overlap of O on the bridge site is higher than that on the top site, with values of 0.79 and 0.76 respectively. This agrees with the earlier conclusion that the bridge site configuration is more stable than the top site configuration for monolayer coverage.

Vibrational frequency analysis of clean and oxygenated diamond (100) surfaces
Surface vibration probes using EELS supports the possibility of top and bridge site co-adsorption. Hossain et al. measured the vibrational energy of the diamond (100) surface with increasing exposure to O atom at 500 K [31]. The C(100)-(2 Â 1) surface has characteristic energy loss peaks at 90 and 152 meV. A loss peak at 215 meV started to appear at low oxygen exposure which has been assigned to carbonyl stretching. At high oxygen exposure (30 L), the loss peaks of clean surface were replaced by loss peaks at 113 meV and 150 meV which have been assigned to carbonyl bending and ether stretching, respectively, while the peak at 215 meV remains. A loss peak at 261 meV also exists which have been assigned to higher-order oxidation states at step edges. The peak at 215 meV is also prominent on a diffuse reflectance infrared Fourier transform (DRIFT) spectra which has been assigned as well to carbonyl stretching. In addition, several peaks between 120 À 180 meV were observed but accurate assignment of these peaks is difficult because of complexity and close proximity to each other [73]. For both studies, the assignment of loss peaks relies mainly on vibration frequency data of small molecules such as 2-adamantanone (C 10 H 14 O) and polycyclic ether dioxabicyclo [3.2.1] octane (C 9 H 16 O 2 ) and did not consider the distinct vibrational features of diamond (100) surface.
To verify the assignment of peaks and to gain theoretical support for top and bridge site co-adsorption, the vibrational frequencies of the clean C(100) surface and C(100) surfaces with 1 ML oxygen coverage on top and bridge sites were calculated. The vibrational eigenvalues were obtained within finite-difference harmonic approximation involving the top-most four carbon layers for clean surface, while the oxygenated surfaces also include the surface O layer. Higher-order oxidation states which gives rise to the peak at 261 meV will not be considered. The results were broadened using Lorentzian with width of 15 meV to match the resolution of EELS spectrometer employed in [31] and is shown in Fig. 12. The EELS spectra for various oxygen exposures from Ref. 31 were redrawn and are shown as well (Fig. 12a). The calculated spectra for the C(100)-(2 Â1) surface have broadbands below $ 150 meV with two main peaks at 80 and 147 meV and a minor peak at 106 meV (Fig. 12b, bottom panel). The two main peaks correspond to the peaks in EELS data. In addition, the calculated spectra also agree well with the high-resolution electron energy loss spectra (HREELS) by Lee and Apai who reported main peaks at 87 and 152 meV and minor peak at 126 meV [74]. Theoretical studies using DFT with local density approximation (LDA) [75] and tight-binding (TB) Hamiltonian [76] reported vibrational spectra with broadbands terminating at slightly higher energies ($ 186 meV) but without the characteristic peaks of 2 Â 1 reconstructed surface of 90 and 152 meV, presumably due to insufficient number of atomic layers in their computational models.
As seen in Fig. 12b (middle and upper panels), simple broadening scheme was not able to reproduce the EELS spectra for diamond surface after oxygen adsorption. In particular, the intensity of the vibration mode corresponding to carbonyl stretching for the C(100)-(1 Â 1):O top surface is very weak (205 and 231 meV). On the other hand, the spectra for the C(100)-(1 Â 1):O bridge surface did not considerably change the features of C(100)-(2 Â1) surface. To address this issue, we calculated the intensity of vibrational modes explicitly. Ab initio calculations have shown that oxygenterminated C(100) surfaces have strong positive electron affinities which gives rise to strong static dipole [35]. This effect has become the basis for the fabrication of various diamond sensors [36,77]. Static dipole could cause intense dynamic dipole when diamond is exposed to time-dependent electric field. Therefore, calculation of vibrational intensities due to dipole scattering mechanism may lead to better vibrational spectra prediction. The intensity data were broadened using Gaussian with FWHM of 8 meV (Fig. 12c). The vibrational modes with the highest intensities are shown in Fig. 13. The consideration of dipole intensities improves the predicted spectra. For the C(100)-(1 Â1):O top surface, peaks can be found at 155 and 231 meV (Fig. 12c, upper panel), which are caused by subsurface phonon bouncing and carbonyl stretching, respectively (Fig. 13a). This affirms the assignment of the EELS peak at 215 meV to carbonyl stretching. On the other hand, for the C(100)-(1 Â1):O bridge surface, peaks can be found at 129 and 153 meV (Fig. 12c, middle panel) which are caused by ether bending and surface phonon bouncing (Fig. 13b). Since bridge site adsorption is the most stable configuration at monolayer, these peaks should correspond to the EELS peaks at 113 and 150 meV observed in high oxygen exposure (Fig. 12a). These results imply that carbonyl and ether can co-exist at high surface coverage. As for the C(100)-(2 Â 1) surface, the peaks found at 122 and 157 meV do not agree with the EEL spectra. This suggests that the EEL spectra of the clean C(100) surface may not be due to the dipole scattering mechanism.
We end this section by proposing a mechanism of oxygen adsorption by synergizing the results of reaction path and vibration spectra calculations with the EELS data. Gas phase O 2 is adsorbed molecularly on the top site and dissociatively on the bridge site. Surface dereconstruction occurs even at low O coverage, and is mainly induced by oxygen adsorption on the top site. This leads to carbonyl formation whose stretching peaks appear at around 215 meV. As the O coverage increases, the isolated epoxide groups transform to ether chain groups and cause surface dereconstruction. This results in the disappearance of peaks at 90 meV and 152 meV and the emergence of peaks at 113 meV and 150 meV. In addition, as the surface reach 1 ML coverage, O atoms of the carbonyl groups along the [011] row move from the top position to the bridge position to form the more stable ether chain. The carbonyl stretching mode did not completely disappear even at high oxygen exposure, suggesting ether chain did not form on some parts of the surface. This could happen along rows with vacancies and other defects.

Oxidation and etching of the diamond (100) surface
The oxidation and etching of the C(100) surface were studied using a larger 2 Â 4 supercell. We have verified that performing simulation on larger supercell will not affect the general conclusions of this section. We started by simulating the desorption of CO and CO 2 from the C(100) with monolayer oxygen coverage on the bridge site and the reaction paths are shown in Fig. 14. Desorption of CO is initiated by cleavage of the CeO bonds of surface ethers and translation of O from the bridge site to the top site to form carbonyl ( Fig. 14A and B). This metastable structure can also be thought of as CO being adsorbed on the bridge site of the second layer atoms (CO bridge ). It is then followed by the near-simultaneous severing of two C-CO bridge bonds giving a total desorption activation energy (E barrier ) of 4:30 eV (Fig. 14BeD). At a distance of 6:9 A perpendicular to the surface, CO is in a physisorbed state (Fig. 14 D) with an energy of 0:14 eV relative to the desorbed state (Fig. 14 E). On the other hand, the desorption of CO 2 proceeds in similar reaction path as the desorption of CO up until the physisorbed state. But just before it fully desorbs, the physisorbed CO bonds with the O atom of the adjacent CO bridge , severs the CeO bond, and desorbs as CO 2 ( Fig. 14F and G). This reaction has an E barrier of 4:88 eV. This extra step makes the likelihood of CO 2 desorption less than CO, agreeing with experiment [31,41]. Thomas et al. reported that the maximum desorption rate for CO and CO 2 both occurs at around 550 C, implying near identical desorption activation energy. This could reflect the similarity of CO and CO 2 desorption process, in agreement with our model [41].
Next, we simulated the succeeding desorption of CO up to the complete etching of the top C layer. The reaction paths and optimized structures are shown in Fig. 15. Table 2 summarizes the desorption activation energies and heat released by the reaction (the energy difference between the transition state and the final   activation energy (reaction A / B and E / F), while succeeding energy barriers are lower (Fig. 15 A-E, E-I). Table 2 shows that reactions A / B and E / F have significantly higher heat released compared to other reactions. Analysis of the reaction path geometries shows that for these two reactions, the initial CO desorption proceeds by the near-simultaneous breaking of two CO bridge -C bond. The point defect left by the initial desorption allows the next CO bridge to break one bond with second layer C and form CO on top (CO top ) structure. The non-simultaneous bond breaking reduces the heat released by the reaction and the desorption activation barrier. This suggests that the point defect functions as nucleation point for CO desorption along [011] direction. The schematic of this process is shown in Fig. 16. This nucleation function was first proposed by John et al. based on their observation that the surface roughness of diamond (100) did not increase even after removal of over 5300 atomic layers, suggesting that the removal rate of rows of atoms is more rapid than the removal rate of layers [32]. Our work is the first theoretical study to support their prediction. Every CO desorbed will leave two dangling bonds on the second layer atoms. But as etching proceeds along [011] direction, these second layer carbon atoms dimerize to attain a more stable surface. A completely etched first layer will leave a second layer with the characteristic 2 Â 1 reconstruction along [011] direction (Fig. 15 I).
The reaction H / I could be viewed as the desorption reaction of isolated CO from the bridge site of C dimer. To gain insights on the bonding mechanism of CO on the surface, we calculated the GPOP of CO molecular orbitals and the corresponding COOP with respect to the bond between CO and the surface. The results are listed on Table 3 while the density of states weighed by GPOP and the COOP curve for frontier orbitals are shown in Fig. 17. The HOMO of gas phase CO is 5s while the LUMO is 2p * , where 2p * is composed of two degenerate orbitals. The calculated gross population of 5s and 2p * are 1:09 and 0:99, respectively. This shows that when CO bonds with diamond (100), an electron from 5s orbital transfers to 2p * orbitals. The DOS weighed by GPOP shows that 5s orbital has no distinct peak which suggests strong hybridization with the bonding states of the C(100) surface. On the other hand, the 2p * orbitals have a sharp occupied peak at À 1:9 eV. The adsorption of CO has total orbital overlap of 0:53 which includes net bonding contributions of 0:27 and 0:32 from 5s and 2p * orbitals, respectively.
The E barrier for the desorption of isolated CO is 1:70 eV. Frenklach et al. modeled CO desorption from the C(100) using a cluster model of C 28 H 48 . In their simulation, the reaction step corresponding to the severing of COeC bond has the highest E barrier of 38:4 kcal=mol (1:66 eV), in reasonable agreement with our result. The activation energy of diamond oxidation was reported in a number of   experiments. For natural diamond and polycrystalline CVD diamond, the typical values were between 2:0 À 2:4 eV [32,39]. Borondoped homoepitaxial C(100) surface oxidized with an activation barrier of 3 eV [31]. Oxidation of single crystal diamond tends to have lower activation energy (1:8 À 1:9 eV) [38,78]. Thomas et al. [41] and Frenklach et al. [33] performed thermal desorption spectroscopy (TDS) and confirmed that the C(100) surface was oxidized by desorption of CO, and reported an activation barrier of 1:91 eV and 1:95 eV, respectively. The most interesting result of TDS data is that the CO desorption peak was very wide, as much as 350 C fullwidth at half maximum (FHWM) for the data of Thomas et al. In light of our simulation results, we propose that the wide peak is due to the existence of multiple peaks corresponding to the wide range of desorption activation energies that we have calculated [41]. The desorption of CO 2 from the oxygenated C(100) surface with defects caused by desorption of first, second, and third CO along the [011] direction was simulated (Fig. 18). The CO 2 desorption following the desorption of the first CO has the highest E barrier of 5:17 eV. It is significantly higher than the E barrier of CO desorption (3:62 eV) from the same surface ( Fig. 15 B/C). The same is true for the two other cases, with E barrier of 2:39 and 3:68 eV, while the E barrier for the corresponding CO desorptions are 1:01 and 1:84 eV, respectively. The results suggest that the primary oxidation mechanism for C(100) is by CO desorption even in the case of defected surfaces.

Stabilization of etched surfaces by oxygen adsorption
A model of oxidized C(100) surface based on methoxyacetone structure has been proposed recently [37]. The surface is characterized by alternating carbonyl and ether, formed on the first and second surface layer, respectively. It is predicted to be more stable than both ether and carbonyl terminated surfaces. The proposed structure could explain the observed co-existence of ether and carbonyl groups on an oxidized C(100) surface. However, the mechanism of its formation has yet to be explained. Analogous structures can be realized when oxygen adsorbs on oxygenterminated diamond surfaces with vacancies formed after CO desorption (Fig. 19). In this section, we simulate the reactions that could lead to the formation of such structures. We begin by simulating the adsorption of gas-phase O 2 on oxygenated surface with two carbon vacancies (Fig. 15C). At a distance of $ 2:1 A from the second layer atoms, a barrier-less ISC occurs as O 2 dissociates and adsorbs on the vacancy sites with E ads ¼ À 7:89 eV, making it more stable than 1 ML ether coverage by an average of À1:32 eV per atom. Next, we simulated the adsorption of O 2 on single-atom deep trough, which is formed by row-wise etching of CO along [011] direction (Fig. 15 E). The reaction path and corresponding geometries for the adsorption of the first and second O 2 are shown in Fig. 20. Unlike the case of O 2 adsorption on two vacant sites, the adsorption on trough is an activated process because the etched surface is already partially stabilized by dimer formation. The calculated reaction path is similar in the case of O 2 adsorption on top of reconstructed surface discussed in Section A. The reaction is preceded by metastable adsorption of triplet O 2 (reaction A / B and D / E), intersystem crossing to singlet O 2 and molecular adsorption (reaction B / C and E / F), and O 2 dissociation and surface dereconstruction (reaction C / D and F / G). Following the adsorption and dissociation of second O 2 , the O atoms move to the bridge site to form the more stable ether chain structure (reaction G / H), giving an average E ads ¼ À5:43 eV. This structure is 0:10 eV more stable than in the case of surface with 1 ML oxygen coverage on bridge site.
Our calculations confirm that surfaces with monolayer oxygen coverage consisting of carbonyl on the first layer and ether on the second layer are structurally stable. Moreover, we demonstrated that these structures can form during oxidation as oxygen is adsorbed on vacant sites left by desorbed species. As a final note, the stabilization of the succeeding layers by oxygen adsorption can explain the stable shallow steps formed during dry oxygen etching of the C(100) surface [38].

Summary of oxidative etching mechanism of diamond (100) surface
In this section, we summarize the results of our calculations and propose a detailed mechanism of oxidative etching of the C(100) surface.  (Fig. 18 B, C, and D). Second, oxygen can adsorb on a vacancy site, bridging two second-layer C atoms and leading to surfaces with carbonyl on the first layer and ether on the second layer. 4. Oxidative etching of the C(100) surface can occur through the desorption of CO and CO 2 . The desorption of CO starts by breaking the CeO bond of surface ether and formation of CO bridge . The desorption of CO 2 follows a similar desorption path but with an extra step of attaching and capturing the O-atom of an adjacent CO bridge , before desorption. The stringent requirement of this extra step will make it less likely to occur compared to direct CO desorption. 5. The desorption of CO leaves a point defect that functions as nucleation point for succeeding desorption, where it allows CO to move from bridge to top position and facilitates a nonsimultaneous bond breaking which lowers the required activation barrier. This made row-wise etching along [011] direction preferable. 6. The CO desorption barrier is around $ 4 eV prior to the creation of the first point defect. Succeeding desorption activation barrier varies between 1 and 3:6 eV, which depends on the resulting surface reconstruction. In the case of isolated CO, the desorption activation barrier is 1:70 eV. The large variation of desorption barriers will result in multiple thermal desorption peaks along a wide range of temperatures. For physical surfaces with non-pristine condition, desorption energies in between the ones we calculated are expected to arise. Hence, the distances between desorption peaks will decrease and the thermal desorption spectra will be wide and continuous with a single peak, similar to the ones observed in the experiments. 7. CO desorption leaves dangling bonds that are strongly stabilized by oxygen adsorption. Etched surfaces along rows in [011] direction will reconstruct and form C-dimers in a similar manner as the reconstruction of the top layer along [011]. Oxygen adsorbs, induces dereconstruction, and forms ether chain on the etched row, with stabilization energies comparable to that of a clean surface. The stabilization of surfaces by reconstruction and oxygen adsorption will prevent the formation of deep etch pits and allow the etching to proceed via shallow steps.

Conclusions
In this study, the oxidation process of the C(100) surface, from the adsorption of gas phase O 2 , to the etching of the first surface layer and its succeeding stabilization, have been investigated using van der Waals-corrected DFT calculations. Initially, triplet O 2 is adsorbed on the reconstructed C(100) surface in metastable state, and at finite temperature, it passes through the triplet-singlet intersystem crossing point several times. Switching to the singlet reaction path leads to a lower energy barrier and more stable final chemisorbed structure, characterized by molecular adsorption on the top site or dissociative adsorption on the bridge site. At low surface coverage, the severing of C-dimer bonds and dereconstruction of the C(100) surface from 2 Â 1 back to 1 Â 1 periodicity is induced by the dissociation of O 2 on top site forming carbonyl structure. At higher surface coverage, surface dereconstruction is further promoted by the formation of the more stable ether chains. For full monolayer oxygen coverage, oxygen prefers to be adsorbed on the bridge site. However, since ether chains will not form next to surface defects, carbonyl species could exist in non-ideal condition. In addition, oxygen strongly adsorbs in vacancies to form ether structure in the second layer. This explains the existence of both ether and carbonyl peaks in EELS experiments. Oxygen coverage beyond 1 ML in form of backbond adsorption is not energetically stable. The oxidative etching of the C(100) surface predominantly occurs through the desorption of CO. Initial desorption of CO requires high activation energy, but with the creation of point defect that functions as nucleation site, successive etching along the preferred [011] direction (perpendicular to the C-dimer) requires lesser activation energy. The etched surface is stabilized through dimerization of subsurface atoms. This leads to a wide-range of activation energies which explains the broad CO thermal desorption peak observed in experiments. Finally, single-atom-deep etched trough is further stabilized by oxygen adsorption. This leads to stable step surfaces and provides theoretical support for a layer-by-layer oxidative etching mechanism.
CRediT authorship contribution statement

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.