Unveiling surface dynamics: in situ oxidation of defective WS₂

Abstract

Applications of transition-metal dichalcogenides (TMDs) are affected by defects and oxidation in air. In this work, we clarify the relationship between oxidation dynamics and O₂ availability for highly defective (and therefore reactive) surfaces of WS₂ crystals. Grazing incidence Ar⁺ sputtering was used to induce a significant concentration of S vacancies in the sample, rendering it highly susceptible to oxidative degradation. In this paper we observe that oxidation occurs slowly under low O₂ pressures (<10⁻⁴ mbar) due to reduced O₂-vacancy interactions. At higher O₂ pressures, the reaction progresses rapidly, as tracked by changes in the oxidation state of W using XPS. The density functional theory calculations support the experimentally observed changes in the oxidation state of W after sputtering and oxidation. They provide the mechanisms of O₂ dissociation on S vacancy clusters, demonstrating that the reaction barrier depends on the coordination of surface W atoms. These results can be useful for protecting samples from degradation in device applications.

---

1 Introduction

The outstanding physical and chemical properties of 2D transition metal dichalcogenides (TMDs) have attracted intense interest in the last decade not least due to their potential applications as sensors, electrocatalysts for the hydrogen evolution reaction (HER), as well as in microelectronics, as reviewed in e.g.^1–3 However, the oxygen reactivity and stability in air of TMD films are still not fully understood, hampering the success of these applications.^4,5 TMD oxidation has been shown to affect their electronic properties,^6–8 induce p-type doping, and affect optical characteristics.⁹ It also leads to the formation of islands on the surface¹⁰ and causes cracking of the material, among other detrimental effects that limit potential applications. Unlike instantaneous damage observed in monolayer black phosphorus,¹¹ oxidation of TMDs occurs gradually over months.¹² To protect TMDs from oxidation, hexagonal boron nitride (hBN) encapsulation is often employed. This significantly enhances their properties,^13–15 highlighting the challenges posed by even small exposure to the atmosphere for device fabrication.

The reactivity of TMDs with oxygen gas has been studied in several publications. Density functional theory (DFT) simulations suggest that the reaction of O₂ with pristine TMD surfaces is highly improbable due to high reaction barriers.^16,17 However, the barrier for the reaction is significantly reduced from 1.59 eV to 0.8 eV in MoS₂ when O₂ interacts with the surface sulphur (S) vacancies.^4,18,19 Scanning Tunnelling Microscopy (STM) experiments reveal a discrepancy with the common assumption that S vacancies are the predominant defects in TMDs; instead, atomic O substitution is observed more frequently after atmospheric exposure.²⁰ Despite these advances, a full understanding of the role of S vacancies in the mechanism of TMD oxidation is still lacking. Further progress can be achieved via controlled creation of surface S vacancies using Ar⁺ sputtering or other methods. For example, in recent work, it was found that oxidation of WS₂ surface occurs more rapidly after controlled Ar⁺ sputtering when the defect concentration (presumed S vacancies) exceeds 10%.²¹ DFT calculations predicted a significantly lower barrier for molecular O₂ dissociation on a S divacancy than on a monovacancy.²¹

In this paper, we further investigate the oxidation dynamics of WS₂ surfaces with high defect densities induced by Ar⁺ sputtering by systematically changing the oxygen pressure and monitoring how controlled oxygen exposure affects surface oxidation using in situ X-ray photoelectron spectroscopy (XPS). Our findings show that, while oxidation remains minimal when O₂ pressure is below 10⁻⁴ mbar, it accelerates rapidly at higher oxygen pressures. Our density functional theory (DFT) calculations show that larger vacancy clusters significantly lower the barrier for O₂ dissociation by reducing the effective charge on W atoms in defect regions. Furthermore, we explore the mechanisms of dissociation of O₂ in singlet and triplet spin states and demonstrate how electron transfer from the surface to the O₂ molecule facilitates dissociation and leads to incorporation of O atoms at various lattice positions.

2 Methods

2.1 Sample preparation

A WS₂ bulk crystal purchased from HQgraphene was used for all experiments. The upper layer of the surface was removed in situ using tape under high vacuum conditions (10⁻⁷ mbar). This ensured a clean and uncontaminated “pristine” surface for exposure. A bulk sample ensures substrate effects such as, charge transfer are limited during the experiment. Subsequently, the sample remained strictly at 10⁻⁹ mbar until the O₂ exposure steps. The sample was then subjected to sputtering with an argon beam, utilising an ion energy of 3 keV and an incident angle of 75° relative to the normal of the sample plane. To confirm uniform sputtering and sample homogeneity, W 4f core level XPS spectra were acquired at various areas of the sample.

O₂ exposure was conducted in the load lock, excluding other gases or atoms present in the laboratory atmosphere, using a liquid nitrogen cold trap. This effectively eliminated the effects of gases such as carbon dioxide, water vapour, and nitrogen. The sample was exposed to O₂ at various intervals at increasing pressures as described in Table 1 resulting in a higher concentration of O₂ molecules per unit of surface area of the sample with each exposure step. The final exposure step was performed at laboratory ambient conditions, in the conditions the sample was exposed to all gases. All measurements and sample preparation steps were performed at room temperature (300 K).

Table 1 O₂ dosing amounts at each exposure step

Exposure step	Total exposure (min)	O2 pressure (mbar)
1	2	5.00 × 10^−6
2	22	5.00 × 10^−6
3	42	5.00 × 10^−6
4	112	5.00 × 10^−6
5	232	5.00 × 10^−5
6	262	5.00 × 10^−4
7	292	5.00 × 10^−4
8	322	1.00 × 10^−2
9	352	1.00 × 10^−1
10	382	1.00 × 10^3
11	412	1.00 × 10^3
12	442	Ambient

---

2.2 XPS

XPS data was acquired using an X-ray photo-electron spectrometer (PREVAC) with a 150 W monochromated Al Kα laboratory source, where the photon energy of the X-rays is 1486.6 eV, with a beam spot size of ∼100 μm and 10⁻⁹ mbar vacuum. The core level spectra were acquired with a pass energy of 100 eV and a 0.1 eV step size, while the survey spectra were collected with a pass energy of 200 eV and 0.5 eV step size. The photoelectrons were collected at normal emission angle, the X-ray incident angle at 60° to the surface normal. All spectra were calibrated to the Au 4f core level of an in situ annealed gold reference sample, where Au 4f_7/2 = 84.0 eV.²²

2.3 XPS data fitting

The core level peaks were fitted using CasaXPS with the GL function for the line shape, incorporating a 70% Gaussian contribution and a Shirley background. Constraints for the W 4f peaks were set as follows: a 2.15 eV spin–orbit splitting, a 1 : 0.75 ratio of area of the W 4f_7/2 to W 4f_5/2 peaks, and a 5.69 eV splitting between the W 4f_7/2 and W 5p_3/2, all based on reference pristine spectra. Additionally, the full width at half maximum (FWHM) of all W 4f peaks was assumed to be the same. The S 2p peaks were fitted assuming a ratio of 1 : 0.5 of S 2p_3/2 and S 2p_1/2 components, and the FWHM of both components was set to be equal.

2.4 DFT calculations

All calculations were carried out using the CP2K code²³ within the Gaussian and plane wave (GPW) formalism²⁴ and the GGA Perdew–Burke–Ernzerhof (PBE) functional.²⁵ A symmetrically extended 6 × 6 × 2 supercell containing 432 atoms was used for the calculations of defect clusters with up to 3 vacancies. For the 7-vacancy cluster, a 9 × 9 × 2 supercell containing 972 atoms was used to minimise strain. All calculations were performed at the Γ point. The symmetric supercell extensions fold important k points into the Γ point and allow us to account for defect-induced lattice distortions.²⁶ Cell parameters were fixed in the defect calculations, with the lattice vectors a and c set to 3.15 Å and 12.08 Å in both the 9 × 9 × 2 and 6 × 6 × 2 supercells. The formation energies for the defects calculated with both supercell sizes differ by less than 0.02 eV.

Periodic 3D boundary conditions with a vacuum gap of 30 Å between slabs were used to ensure that there was no interaction between the slabs. We used the MOLOPT-DZVP-SR basis set²⁷ and GTH pseudopotentials²⁸ optimised for PBE. The van der Waals interaction was taken into account using a DFT-D3 dispersion correction by Grimme et al.²⁹ The cutoff and relative cutoff used in all calculations were 600 and 60 Ry, respectively. All effective charges are in units of |e| and calculated using Bader population analysis software by the Henkelman group.³⁰ Kohn–Sham (KS) one-electron energy differences were used to estimate band gaps. Climbing image nudged elastic band calculations (CI-NEB)³¹ with 7 replicas were used to calculate the barriers for the dissociation of the O₂ molecule and migration of the S vacancy. AIMD calculations were performed at 300 K with a time step of 2 fs in the NVT ensemble using a CSVR thermostat.³²

3 Results

3.1 Characterisation of WS₂

To establish a baseline measurement of the pristine surface of WS₂, the sample was cleaved in situ in the load lock at 10⁻⁷ mbar, using Scotch tape, as depicted in Fig. 1 and detailed in section 2. The cleaved sample was then transferred into the UHV XPS system for measurements. Fig. 1 shows a schematic of the overall sample treatment in this work and the resulting W 4f spectral changes, these will be discussed in more detail below. The survey spectra of the sample (Fig. S1†) demonstrate that the WS₂ surface is clean, with no detectable oxygen and only trace amounts of carbon. The W 4f_7/2 peak is observed at 33.2 eV, which corresponds to W in the nominal +4 oxidation state in WS₂, as illustrated in Fig. 1. The peak exhibits a full width at half maximum (FWHM) of 0.6 eV and a spin–orbit splitting (SOS) of 2.15 eV, consistent with previously reported values.³³ The S 2p_3/2 peak is located at 162.9 eV, consistent with S in TMD compounds. The FWHM of this peak is 0.7 eV, with a SOS of the S 2p_3/2 and S 2p_1/2 levels being 1.2 eV in agreement with the results of ref. 34.

Fig. 1 Overall process of sample treatment during this work with corresponding W 4f core level spectra. The graphic illustrates the treatment of the sample. Initially, bulk WS₂ is cleaved in situ. Subsequently, the sample undergoes Ar⁺ sputtering to induce defects such as S vacancies (V_S), after which the sample undergoes oxidation. The W 4f core level spectra are then measured and shown below for each sample treatment: pristine, sputtered, and oxidised WS₂. Each W 4f core level corresponding to different W environments is fitted with a doublet for the 4f_7/2 and 4f_5/2 peaks for the pristine WS₂ contribution, the shoulder arising from sputtering (WS_2−x) and the oxide contribution (WO₃), as well as the 5p_3/2 peak.

After in situ exfoliation, the sample was exposed to laboratory atmosphere for 80 hours and subsequently re-measured. XPS analysis of core levels revealed no discernible changes in the chemical states of W and S (Fig. S2†), indicating relatively high stability of the pristine surface with only a slight increase in the relative signal intensity of C 1s and O 1s core levels post-exposure from adsorption of contaminants on the surface. However, defective WS₂ degrades relatively rapidly according to ref. 8 and 35. In particular, Bussolotti et al.²¹ have recently demonstrated that samples with a defect density of 7% or below exhibit minimal oxidation for several hours in a laboratory atmosphere. However, samples with a defect density exceeding 20% undergo significant oxidation. Consequently, in the presence of a high defect density, monitoring of the oxidation of the sample becomes feasible due to pronounced changes in the electronic structure. In this work, the sample was sputtered until a 20% reduction of the S 2p core spectrum area was observed, which allows us to systematically track surface oxidation.

Sputtering results in four spectral changes: (i) a BE shift of all levels; (ii) appearance of a shoulder in the W 4f core-level at 31.8 eV; (iii) a decrease in the area of the S 2p core-level; (iv) and an increase of the FWHM of the W 4f and S 2p core levels. These effects are illustrated in Fig. 2, which depicts the W 4f and S 2p core level spectra before and after sputtering, as shown in panels (a) and (c), respectively. A shift of 0.6 eV in BE for all levels is the result of the shift in the Fermi level of the sample, which all levels are measured relative to in XPS spectra. The peaks corresponding to WS₂ shift to 32.6 eV and 162.3 eV for W 4f_7/2 and S 2p_3/2, respectively. A shoulder that appears in the W 4f core-level spectra at 31.8 eV signifies a change in the W environment, such as the removal of S atoms. We refer to the origin of this shoulder caused by sputtering as a WS_2−x species, because the exact nature of the surface is not yet known.²¹ Peaks at lower binding energies are characteristic of the element undergoing reduction, attributed to the enhanced electron density surrounding the element and alterations in the Coulomb interactions between the photoemitted electron and ion core.

Fig. 2 Changes in core level spectra during various sample treatment stages. (a) The W 4f core level before and after sputtering showing the emergence of a shoulder at a lower BE, broadening of the line shape (from 0.6 to 1 eV) and a 0.6 eV shift in BE. (b) The ratio of total area of the W 4f core level before (P which stands for Pristine) and after sputtering (S which stands for Sputtered) and during the different O₂ exposures, relative to the pristine sample peak area of WS₂. The total area of W remains relatively constant around ±0.05 (10%). (c) The S 2p core level before and after sputtering showing the broadening of the line shape (from 0.7 to 1 eV) and a 0.6 eV shift in BE. (d) The ratio of total area of the S 2p core level before and after sputtering and during the different O₂ exposures, relative to the pristine sample peak area of WS₂. The area decreases by 0.2/0.25 (20/25%) and then by a further 10% in the final preparation step. Plots of measured XPS core levels during different O₂ exposures. (e) W 4f (f) O 1s (g) S 2p and (h) C 1s. All spectra have been normalised to the total relative area of W 4f.

The summary of the total area of the core level spectra relative to the pristine spectrum is shown in Fig. 2(b) and (d) for W 4f and S 2p, respectively. The first two points illustrate the changes in pristine and sputtered samples. It can be seen that, while the total area of the W peak remains relatively constant (within ±10% as shown in Fig. S3(a)†), the area of the S core peak decreases by 20/25% after sputtering (Fig. S3(b)†). We hypothesise that defects induced by sputtering are predominantly S vacancies. This is supported by the reduction in the S 2p area demonstrated in ref. 36. for sputtering under the same conditions. This preferential sputtering of light elements, which is more prominent at incident angles around 75°, is also seen in other oxides, such as Ta₂O₅,³⁷ whilst sputtering of W metal at this angle has very low sputtering yields,^38,39 as has also been shown theoretically in TMDs.^40–42

According to ref. 43 the S-vacancy concentration (C) on the sample surface can be expressed as:

where A and A′ are the intensity of the W 4f peak of the stoichiometric WS₂ and substoichiometric WS_2−x, respectively, and the pre-factor 1/3 reflects the atomic coordination around the vacancy. Based on ref. 41, the S vacancies and the related WS_2−x phase is mainly located within the surface layer. In contrast, the W 4f signal of the WS₂ phase also includes signals from the inner WS₂ layer due to the finite electron mean free path of electron (∼2 nm) at the investigated photoelectron kinetic energy values (∼1400 eV (ref. 44)). As such, this equation provides a lower limit of the surface defect concentration. Based on the measured XPS intensity in the W 4f spectra of the sputtered sample a value of C > 8% is obtained.

Moreover, the W 4f and S 2p peaks in WS₂ also exhibit an increase in FWHM to 1 eV, as determined by fitting. The final line shape results from a combination of factors including broadening due to core hole lifetime (Lorentzian contribution), spectrometer broadening, charging effects, and geometric effects (such as disorder due to phonons which is the Gaussian contribution to the lineshape). The core hole lifetime remains constant for a core level of an element in a specific oxidation state, and the spectral broadening also remains constant throughout the measurement. Charging effects that are often observed alongside a shift in the BE of the core level were not present in this work, therefore, this effect can also be excluded. The observed broadening is probably caused by changes in the local structure of the sample and suggests an increase in the disorder in the surface structure, as has been suggested in ref. 40. In particular, no discernible new peaks emerge in the S 2p core level spectra, indicating no significant alteration in the environment surrounding the S atoms.

In the next section we present the results of exposure of the sample after sputtering to O₂ at various time intervals and pressures using the approach described in detail in section 2.

3.2 Oxidation of sputtered WS₂

Initially, no detectable chemical changes occur on the sample surface. However, when the partial pressure of O₂ is increased to more than 10⁻⁴ mbar, the shoulder that was created by sputtering begins to disappear in the W 4f core level, while a new shoulder emerges at 35.8 eV (Fig. 2(e)). This higher BE peak suggests the presence of a W–O bond,⁴⁵ indicating the formation of oxide on the surface. The disappearance of the shoulder at 31.8 eV suggests the passivation by exposure to oxygen of S defects created by sputtering. In the final exposure step, the sample is exposed to atmospheric conditions, resulting in nearly complete oxidation. The sample was exposed to all gases, including water vapour in atmosphere, during the final exposure step. This may affect the oxidation process. The peak at 35.8 eV corresponds to W 4f_7/2 in the +6 oxidation state in WO₃.

Gradual exposure leads to the formation of a thin O₂ and carbon layer on the surface of the sample, which was absent immediately after sputtering, as observed in the O 1s and C 1s core level spectra (Fig. 2(f) and (h), respectively). The emergence of surface oxide is accompanied by a significant increase in the O 1s intensity. The peak is asymmetric, which confirms the presence of at least two different O environments. There is a change in relative peak intensity during exposure, with most of the increase in intensity of the area coming from the lower BE feature of the peak. This will be discussed further below. Meanwhile, during the oxidation process, the S 2p core level remains unchanged (Fig. 2(g)).

At all values of O₂ pressure, the line shape of the C 1s core level remains unchanged. While O 1s exhibits a significant increase in the second (lower BE) peak at 530.8 eV. The emergence of this second peak in the O 1s core level coincides with the appearance of a shoulder in the W 4f core level at higher BE. This suggests that, in addition to the physisorbed O, the presence of which we attribute to the peak at 532.3 eV, O₂ becomes chemically integrated into the WS₂ lattice. This integration further implies the formation of W–O bonds⁴⁶ and the subsequent development of a surface oxide layer. Physisorbed species can be attributed to hydroxide/water from the walls of the UHV chamber and physisorbed O₂.^47,48 The changes for subsequent exposure steps in the individual W 4f core level are illustrated in Fig. 3.

Fig. 3 Examples of deconvoluted W 4f spectra and their relative contributions. (a) The fitted W 4f core level spectra, with different chemical environments labelled including WS₂, WS_2−x, and WO₃. The numbers in the top LHS corner correspond to the preparation step at which the spectrum was collected. (b) Breakdown of the different chemical environments of the W 4f spectra by area, relative to the total area of W 4f_7/2 only, demonstrating the decrease of the sputtered (WS_2−x) peak as the WO₃ peak increases.

The W 4f core levels and some of the corresponding fits for all exposures are shown in Fig. 3(a) (for all fits of the W 4f see Fig. S4†). Initially, the introduction of O₂ induces minor changes in the W spectra (exposure steps 1 to 7). The ratio of the sputtered peak area (W 4f_7/2 in WS_2−x) to the total W 4f area (of the doublet) remains relatively constant (between 0.25 and 0.21) until the eighth step of oxidation. This makes it challenging to discern if there is an overall change in intensity, as summarised in Fig. 3(b). By the 8^th exposure step (pressure of O₂ above 10⁻² mbar), a noticeable decrease in the WS_2−x shoulder is observed. This suggests that W atoms, which have been reduced by sputtering, are reacting with O₂ as their oxidation state increases (S vacancies are being passivated). An oxide layer begins to form and this oxide layer is fitted with another doublet, as described in section 2. The W 5p_3/2 for the oxide peak is omitted due to its low or ‘negligible’ intensity. Subsequently, after the 8^th exposure step, the WS_2−x peak ratio decreases significantly, and the oxide peak ratio becomes significant. By dosing steps 10 and 11, the oxide peak area surpasses the WS_2−x peak area. After the final exposure, the WS_2−x peak area is nearly zero, and the oxide peak has a considerable intensity, indicating the oxidation of most of the reduced W atoms on the surface.

The smaller oxide peak area in the final exposure step, compared to the WS_2−x peak pre-oxidation, indicates the possibility of incomplete conversion of W to the +6 oxidation state. This can occur as a result of the formation of a tungsten oxysulfide species with W in the +4 oxidation state, consistent with previous findings.⁴⁹ The peak for this species is difficult to fit due to the overlap with WS₂ peak and therefore, was omitted in our fitting. The corresponding residual plots for the fits are shown in Fig. S4.† Furthermore, this is accompanied by the rise of the peak at lower BE, 530.8 eV, attributed to lattice O in WO₃.^33,50 The corresponding fitting of the O 1s core level can be seen in Fig. S5.† The rapid increase of the peak at low BE with respect to the peak at high BE is observed, suggesting that the concentration of physisorbed O₂ levels off at higher O₂ doses (as the surface becomes saturated), and the increase in the area of the O 1s core level is mostly due to the formation of surface oxide. We cannot completely exclude the presence of an oxide layer before the 8^th exposure step, but if it is present, it is below the detection limit of our experiment.

A summary of all changes at the core level during oxidation is illustrated in Fig. 4. The disappearance of the WS_2−x peak (Fig. 4(a)) is modelled with two linear fits between exposure steps 0 and 7 (grey line – segment 1) and between 8 and 12 (pink line – segment 2). This approach is justified because the decrease in segment 1 is notably smaller compared to segment 2. The decrease rates obtained from the fit gradient are −0.003 and −0.02 for segments 1 and 2, respectively. This indicates that the decrease in the WS_2−x shoulder is an order of magnitude faster at O₂ pressures of 10⁻² mbar and above. In segment 1, the data has a scatter of ±0.01 due to measurement and fitting uncertainties. Segment 1 exhibits significantly higher scatter, as evidenced by the R² value of 0.47 for the linear fit presented in the caption of the figure. Therefore, within the error of the experiment, the oxidation of defective WS₂ below these O₂ pressures is insignificant. The S 2p core level remains unchanged, with intensity fluctuations ranging from 0.38 to 0.46 (Fig. 4(b)) and error bar overlap. A linear fit yields a gradient of −0.002 and an R² value of 0.09, indicating a weak correlation between the area of the core level and the amount of O₂ to which the sample is exposed. Changes in the area of O 1s (Fig. 4(c)) increase with the same rate as the decrease in WS_2−x when the surface oxide is formed, whilst the C 1s increases rapidly and levels off after the 4^th exposure step (Fig. 4(d)) suggesting no/low lattice incorporation of C.

Fig. 4 Plot of areas of core level peaks at different stages of oxidation normalised to the total W 4f core level area. (a) W 4f_7/2 fitted between the 0 and 7^th, and 8^th and 12^th exposure steps due to the distinct changes in the two regions. The fit parameters are as follows: Segment 1 (0–7, grey line) Gradient: −0.003, R²: 0.47, and Segment 2 (8–12, pink line): Gradient: −0.02, R²: 0.95. (b) S 2p – the linear fit has a slope of −0.002 and R² of 0.09. (c) O 1s fitted between the 0 and 7^th, and 8^th and 12^th exposure steps due to the distinct changes in the two regions. The fit parameters are as follows: Segment 1 (0–7, grey line) Gradient: 0.004, R²: 0.65, and Segment 2 (8–12, pink line): Gradient: 0.02, R²: 0.90. (d) C 1s – fitted with an exponential function y = a·exp(−bx) + c, where a = −0.03, b = 0.49, c = 0.03. The gradients describe the rate at which each peak increases (relative to W) with each exposure step. High R² values are indicative of smaller deviations of the data from the linear fits.

The reaction of O₂ with sputtering-induced S vacancy defects (as described in section 3.1) leads to a reduction in the area of the WS_2−x peak in the W 4f core level, accompanied by a concurrent increase in both the oxide peak of W 4f and the O 1s peak at O₂ pressures of 10⁻² mbar and higher. Previous studies have shown that the barrier for O₂ dissociation at S vacancy dimers is significantly reduced,²¹ suggesting that oxidation is more likely to occur at vacancy clusters (where N = 2 or more). We expect that vacancy clusters of various sizes will be present on the surface post-sputtering.

We note that vacancy clusters and pits have been shown to form in MoS₂ as a result of Ar⁺ sputtering in ref. 51–54. In the following, we use DFT to demonstrate how the effective charge on W varies with the vacancy cluster size, linking this to changes in the oxidation states (and therefore the core level shifts observed in our spectra). Variations in effective charges have previously been directly linked with changes in oxidation states.^55–57 Furthermore, we show that the barrier for O₂ dissociation in vacancy clusters reduces as the cluster size increases and consider the effect of the triplet to singlet transition during the O₂ dissociation.

3.3 Electronic properties of the WS₂ surface and its interaction with O₂

3.3.1 Properties of vacancy clusters. The structures and properties of various S vacancy clusters at the WS₂ surface, ranging from single to seven S vacancies were calculated. They include dimers and trimers (in line, angled, and triangular arrangements), as well as larger clusters shown in Fig. 5(a). In this work we assume the creation of predominantly S vacancies, although we cannot exclude the formation of W vacancies. However, due to differences in the sputtering yield and therefore the increased likelihood of ejecting S atoms during sputtering, as discussed above, our DFT calculations mainly consider S vacancy clusters.

Fig. 5 S vacancy complexes calculated in this work and their in text notation. (a) Models of calculated vacancies and (b) corresponding effective charges (Q_eff) on surrounding W atoms. The charge on S atoms remains relatively unchanged upon vacancy creation. This corresponds to the additional shoulder at low BE in the W core level spectra, while the S environment remains unaffected. For the pristine structure, the reference W Q_eff is 1.24. Images have been cropped for clarity and do not represent the whole system used in calculations. Up to 3 S vacancies have been calculated using a 6 × 6 × 2 supercell and 9 × 9 × 2 supercell for V7 cluster.

Using DFT and Bader charge analysis, we determined the effective charges of W atoms embedded into different S vacancy clusters. The change in the effective charge of the W atoms is proportional to the number of missing bonds and the size of the vacancy cluster. These charges are 1.04, 0.93, 0.73|e| when 1, 2 and 3 S vacancies surround a W atom, respectively (Fig. 5(b)). The effective charge on W in the middle of the 7V cluster decreases further to 0.65|e|, suggesting that the W atoms inside the clusters may exhibit increased reactivity.

The observed reduction in the oxidation state leads to a negative shift of the core level and can contribute to shoulder formation observed at 31.8 eV in the W 4f spectra. Larger vacancy clusters result in more empty in-gap states (Fig. S6†). Consequently, electrons from the valence band are more likely to be excited into the vacancy-induced empty states, which could explain the 0.6 eV shift in the binding energy (BE) observed post-sputtering. The shift in XPS spectra is attributed to a change in the Fermi-level position of WS₂.

However, cluster formation is not accompanied by an energy gain. The energy change per vacancy is calculated as , where E^f_N is the formation energy for the N isolated vacancies and E^f_complex that of the complex. For 2 S vacancies (2V), 3 S vacancies (line – L3V –, V-shaped – V3V –, and triangle – T3V – configurations), and the 7 S vacancy complex (7V), we obtain E_bind values of −0.01, 0.00, −0.02, −0.12, and −0.24 eV, respectively, as shown in Table S1.† The cell parameters used for the calculations can be found in section 2.4. Negative values indicate that binding is unfavourable. For clusters with more than three S vacancies, other cluster shapes are possible, which may have varying E_bind values. However, we expect that the thermodynamic stability of the clusters will not play a significant role in our conclusions because Ar⁺ sputtering results in the formation of clusters of varying and random configurations, which are unlikely to rearrange, as explained below.

3.3.2 Sulphur vacancy rearrangement. To better understand the surface morphology and vacancy distribution after sputtering, we investigated the barriers of S vacancy migration. The various processes considered are illustrated in Fig. 6(a), with their corresponding barriers shown in Fig. 6(b). Most vacancies are expected to reside in the upper S-plane of the top layer due to the high incidence angle used during sputtering.⁴⁰ However, if any vacancies form in the lower S-plane of the top layer, their diffusion into the upper layer has low probability under experimental conditions (300 K) due to the high barrier of 4.93 eV. The diffusion of single vacancies across the surface, as well as dimer rearrangement processes, have barriers of 3.01 eV and 2.54 eV, respectively, and will also proceed slowly.

Fig. 6 S vacancy diffusion barriers on the WS₂ surface. (a) Diffusion trajectories in which a S atoms migrates from the position represented by the green ball to the position represented by the yellow ball, the S vacancy sites are represented by the dotted circle and (b) corresponding barriers for various processes.

It should be noted that these barriers are higher than previously reported for monolayer MoS₂, consistent with previous work.⁵⁸ The difference arises because the bond dissociation energy (BDE) of WS₂ is higher than that of MoS₂ by around 1 eV, with a W–S and Mo–S BDE of 4.94⁵⁹ and 3.94⁶⁰ eV, respectively, as measured by resonance two-photon spectroscopy. The migration barriers for S vacancies on the WS₂ surface are relatively high, because strong W–S bonds must be broken in order for S atoms to migrate from one site to another. The interlayer diffusion barrier (process 1 in Fig. 6(a)) is nearly 4.93 eV, which can be attributed to the bond breaking without any compensatory bond formation in the transition state. Therefore, it is unlikely that the vacancies will migrate across the surface at room temperature within the timescale relevant to our experiments. Furthermore, there is no thermodynamic drive for vacancy clustering, as described above.

However, at high O₂ pressure, rapid oxidation of the sputtered samples is observed experimentally, suggesting that Ar⁺ sputtering creates reaction sites. Previous research has demonstrated that O₂ dissociates with a relatively low barrier upon direct collision with the S vacancy dimer.²¹ Below, we investigate this trend further by studying reaction of O₂ molecule with an S vacancy dimer, trimer and then heptamer.

3.3.3 O₂ molecule reaction with S vacancy dimer and trimer. On the pristine surface, O₂ physisorbs in a triplet state, with a binding energy of 0.14 eV. This is similar to previously reported values for TMDs^19,61–63 which are all broadly within 0.05–0.23 eV. An ab initio molecular dynamics (AIMD) simulation of the surface with O₂ molecules shows that at 300 K these molecules readily desorb and diffuse around the surface (Video S1†). The dissociation barrier of O₂ on a pristine WS₂ surface has been calculated at 1.76 eV (ref. 4) and at a single S vacancy at 0.86 eV.¹⁹ We find that the dissociation of O₂ is affected by the spin state of O₂ molecule and the effective charge on W in the S vacancy clusters shown in Fig. 5(b). For example, the dissociation barriers are 0.17 and 0.29 eV for a 2V in a singlet and a triplet spin state, respectively. However, the probability that the O₂ molecule will encounter 2V at the correct geometry for the dissociation to occur is low and should result in a low cross section of this reaction. As O₂ molecules easily diffuse around the surface, it is more likely that dissociation will occur at larger vacancy clusters.

Next, we consider the mechanism of O₂ dissociation on a L3V cluster shown in Fig. 5(a). This reaction also does not occur spontaneously. The O₂ adsorption energy on L3V is 0.15 eV. The geometry of O₂ during adsorption is the same in the singlet and triplet states, with O₂ adsorbing at an 15° angle to the surface plane (nearly parallel), angled toward the vacancy cluster, which is similar to previously reported ones (see ref. 64). Because of weak physisorption, the molecule is not ‘anchored’ to the surface and likely moves around, occasionally desorbing. The triplet spin state is the most stable adsorption configuration of O₂, with an energy 0.34 eV lower than that of the most stable singlet state, which is typical for O₂ molecular physisorption. We note that the singlet–triplet splitting is underestimated in our calculations by 0.6 eV with respect to the experimentally reported values.⁶⁵ In the triplet state, the Highest Occupied Molecular Orbital (HOMO) of the WS₂ surface and the Lowest Unoccupied Molecular Orbital (LUMO) of O₂ in the β spin state are 0.46 eV apart. This adsorption is characterised by the interaction of the W d states and O p state as reported for metallic surfaces.⁶⁶ In the singlet state, the HOMO–LUMO splitting is reduced to 0.12 eV, as seen in the projected density of states (pDOS) in Fig. 7(a).

Fig. 7 The process of O₂ dissociation in singlet and triplet states. (a) Comparison of the electronic structures of physisorbed O₂ on the surface in triplet and singlet spin states. The HOMO (surface states) to LUMO (π* oxygen orbital) gap is 0.46 eV for the triplet state and 0.12 eV for the singlet state. (b) Changes in the dissociation barrier for O₂ into the lattice for both spin states, with a lower barrier observed for dissociation of singlet O₂ into singlet O incorporated into the lattice. All energies are referenced relative to triplet O₂ physisorbed on the surface.

The energy barriers for O₂ dissociation on L3V in the singlet and triplet spin states are 0.19 eV and 0.33 eV, respectively (Fig. 7(b)). The final dissociated configuration is by 0.85 eV more stable in the singlet than in the triplet state. The lower barrier in the singlet state is expected because the occupied surface states are closer to the unoccupied LUMO of O₂ (π* orbital), facilitating dissociation, driven by electron transfer from the surface to the molecule, similar to the Ag(110) surface.^67,68 A lower dissociation barrier has also been reported on the Ni(111) surface for the singlet spin state.⁶⁹ The reaction is highly exothermic with the energy of the most stable final geometry, where both O atoms are embedded in the L3V and bonded to W in the lattice, being by 8.62 eV lower than when O₂ is in the gas phase. This large energy gain can be understood if one takes into account that the O atom incorporation energy into S vacancy is 4.42 eV. Similar barriers and reaction energies of O₂ with S divacancies have been predicted in MoS₂.¹⁷ The large amount of energy released will dissipate into generating phonons, electron–hole pairs (see also discussion in ref. 66) as well as diffusion and desorption of atoms and molecules, particularly SO₂.¹⁷

We note that O₂ dissociation on V3V (see Fig. 5(a)) is similar to that of L3V (and 2V) with a barrier of 0.13 and 0.30 eV in singlet and triplet states, respectively. The dissociation on T3V in singlet and triplet spin states is spontaneous with no barrier. We relate this to a much lower effective charge on the exposed W in this cluster, as can be seen in Fig. 5(b). This shows that the exposure of W atoms and their electron population can be important factors for O₂ dissociation.

To attain a deeper insight into the O₂ dissociation mechanism, in Fig. 8, we analyse the properties of each replica from the NEB calculation of O₂ dissociation on L3V in both singlet and triplet spin states. Fig. 8(a) shows the energy of the system as a function of the distance between the O₂ molecule and the surface, defined as the distance from the centre of mass of the O₂ to the plane of W atoms in the upper layer. Initially, the O₂ molecule is physisorbed above the vacancy at a distance of 3.83 Å and 3.70 Å from the W atoms in the singlet and triplet states, respectively. The length of the bond of O₂ is 1.23 Å in both spin states, as shown in Fig. 8(b).

Fig. 8 The relative energy as a function of system properties in triplet and singlet spin states during O₂ dissociation. (a) The distance of the centre of mass (COM) of the O₂ molecule from the surface, defined by the plane of W atoms in the upper layer. The final distance for O incorporated into the lattice is 1.00 and 0.96 Å in the singlet and triplet spin state, respectively. (b) The distance between the two O atoms in the O₂ molecule. The final values are 3.09 Å for the singlet state and 3.06 Å for the triplet state after dissociation. (c) The effective charge on both O atoms. A total of 2 electrons is transferred from the surface to the O atoms after incorporation into the lattice. The dotted line indicates the properties of the system at which −0.76|e| and −0.71|e| is transferred to O₂ in the singlet and triplet spin state, respectively. Corresponding to an O₂⁻ intermediate.

In the triplet state, when the O₂ molecule reaches 1.98 Å from the W plane within the vacancy cluster (dotted line in Fig. 8(a)), 0.71|e| of charge is transferred, as seen in Fig. 8(c) (dotted line). This transfer elongates the O₂ bond to 1.43 Å consistent with the O₂⁻ species, which typically has a bond length of 1.32 Å when chemisorbed on a Pt(111) surface.⁷⁰ Previous studies suggest that metal oxidation processes proceed through an O₂⁻ intermediate.^70,71

In the dissociated configuration, where O atoms become embedded in the lattice, a total of 2|e| are transferred to O atoms. Once fully dissociated, the bond distance of the W–O bond is 2.08 Å. After dissociation, the O atoms are separated by about 3.1 Å.

This suggests the following reaction mechanism: the O₂ molecule, in the triplet state, randomly encounters a vacancy cluster. As the molecule approaches the surface, at 1.98 Å electron transfer from the surface to the O₂ molecule begins, forming a bond between the molecule and the W atoms. The occupation of the anti-bonding orbital π* weakens the O–O bond, increasing its length to 1.43 Å. The bond then weakens further, leading to dissociation and occupation of the S sites driven by the exothermic nature of the reaction.

The low barrier for the O₂ dissociation on 3V clusters implies that increasing the number and size of vacancy clusters, as well as the O₂ pressure, will enhance the probability of oxidation, consistent with experimental observations of higher oxide formation rates at increased oxygen pressure. However, the high mobility of O₂ on surface terraces suggests that bigger vacancy clusters could provide more efficient dissociation sites because of the larger number of exposed W atoms with lower effective charges. In the following, we consider an example of an S heptavacancy (7V) cluster exposing 7 W atoms.

3.3.4 O₂ molecule reaction with S vacancy heptamer. To explore the various configurations of O atoms inside 7V we used a 9 × 9 × 2 supercell. These clusters can adopt numerous configurations, but evaluating their relative stabilities is beyond the scope of this work and we used a symmetric configuration shown in Fig. 5(a). As can be seen in Fig. 5(b), the effective charges of the exposed W atoms are further reduced with respect to the L3V and the perfect lattice. We observe that there is spontaneous dissociation in some configurations of O₂ when it approaches the cluster about 3.66 Å away from the W atom plane (in the z direction). At least one O atom must interact with the high electron density W atoms of the vacancy cluster for dissociation to occur spontaneously during the geometry optimisation. As in the case of smaller clusters, the dissociation is accompanied by the electron transfer from the WS₂ layer.

When O₂ reacts with the vacancy clusters, during the process of the charge transfer from the W atoms, O atoms passivate in-gap states (Fig. S7(a)†). This increases the effective charge of W atoms (Fig. S7(b)†), lowering their electron density compared to the pristine structure (1.37|e| for one O substitution, up to 1.66|e| for 7 O substitutions, versus 1.24|e| in the pristine case). This change, driven by O's higher electronegativity, aligns with the appearance of the oxide peak at 35.8 eV during oxidation.

Fig. 9(a) illustrates the possible configurations of oxygen atoms after O₂ dissociation within the cluster. We consider only configurations where the O atoms are adjacent to each other, with a maximum separation of 3.1 Å (which is approximately the distance between two lattice positions of S), corresponding to their equilibrium distance for incorporation into the lattice during O₂ dissociation. The positions analysed include lattice sites (filled red circles), interstitial positions (dotted empty circles), and one configuration in which an oxygen atom occupies a position inside the cluster while another sits atop a S atom (green circle). Positions marked with red crosses were excluded, as oxygen atoms in these positions either caused significant displacement in the surrounding lattice due to proximity to S atoms or migrated into the more stable position atop S (green circle). Although we cannot entirely rule out such processes in experiments, they are likely less favourable.

Fig. 9 Configurations of O atoms in a 7 S vacancy cluster. Oxygen positions are shown in relation to lattice sites (filled red circles), interstitial sites (dotted empty circles), and positions on top of S atoms (green circle). Red crosses mark excluded positions due to significant lattice distortion or steric hindrance. Symmetrically equivalent configurations are denoted by an equals sign.

Fig. 9(b) summarises the relative energies of the various oxygen positions compared to the most stable configuration at the position 5 in the centre of the cluster in connection with other lattice positions (1, 2, 4, 6, 9, 10). Lattice sites are inherently more stable as a result of reduced steric hindrance. The next most stable configurations involve O atoms occupying outer lattice positions (e.g., O1 & O2 in positions 1 & 2), which are 0.2 eV higher in energy due to the proximity of both oxygen atoms to S atoms in the lattice. Configurations where one oxygen atom occupies a lattice site and the other occupies an interstitial site (e.g., positions 1 & 3) are 2.3 eV higher in energy. When both oxygen atoms occupy interstitial sites (e.g., positions 3 & 7), the configuration is 2.4 eV higher. Finally, if O₂ interacts with the edge of the vacancy cluster, it can dissociate into a lattice position and one site atop S (e.g., positions 9 & 10), with this configuration having an energy of 4.6 eV higher, corresponding to a single substitution of the S vacancy. We note that, despite the difference of 4.6 eV between the relative stabilities of the oxygen atoms dissociated in the cluster, all configurations are more stable than molecular O₂ physisorbed on the surface.

4 Discussion and conclusions

We used grazing incident Ar⁺ sputtering to accelerate oxidation of WS₂ layers. Our findings reveal that defects induced by sputtering serve as highly reactive sites for oxidation and the reaction progresses rapidly under moderate O₂ pressures (>10⁻⁴ mbar). This behaviour marks a departure from previous studies, which primarily focused on pristine or less damaged samples where oxidation occurs gradually over extended periods. Here, we demonstrate that high defect densities significantly alter oxidation kinetics and that the oxidation process can be slowed at O₂ pressures lower than 10⁻⁴ mbar. We note that this threshold is not a universal constant and can vary depending on factors such as the density of S vacancies, surface roughness, sample temperature, and experimental conditions. While further work is needed to explore variability across different samples and setups, this threshold highlights the critical role of oxygen availability in defect-driven oxidation.

XPS measurements revealed that post-sputtering the W 4f core levels exhibited a 0.6 eV shift, reflecting a shift in the Fermi level due to overlapping defect states and changes in the surface morphology or changes in electronic structure as a result of adsorbate removal. The presence of a shoulder at 31.8 eV in the W 4f spectra is consistent with the formation of reduced W atoms, indicative of S removal. A 20–25% reduction in the S 2p peak area after sputtering, indicates the formation of a high density of S vacancies. Our investigation suggests that the sputtered surface is likely to possess S vacancy clusters of different sizes. Upon O₂ exposure, this shoulder decreased significantly, while a peak at 35.8 eV, associated with W in the +6 oxidation state (WO₃), emerged. This transformation was more pronounced at O₂ pressures greater than 10⁻² mbar, where nearly complete surface oxidation was observed, with the WS_2−x peak area decreasing by approximately 0.02 normalised area units per exposure step as opposed to 0.003 in the low O₂ pressure regime.

Valerius et al. have shown that sputtering with Xe⁺ ions at a grazing incidence angle (same as in this work) primarily removes top S atoms, causing reversible disorder.⁴⁰ This is consistent with the broadening of the W 4f and S 2p core level peaks observed in our study, indicating surface disorder. Prior studies have shown the preferential removal of S atoms during sputtering in TMDs due to their lower atomic weight,^72,73 as well as the low sputtering yield of W metal at a 75° incident angle of Ar⁺ ions.^38,39 The Fermi level shift after sputtering can occur as a result of many overlapping defect states from large S vacancy clusters.⁷⁴ In this work we assume that the vacancies generated under these sputtering conditions are primarily located in the uppermost layer of the sample. However, we cannot entirely rule out the presence of vacancies in deeper layers. Consequently, we also assumed that any subsequent oxidation primarily occurs in the top layer. Further investigation is required to determine the precise surface morphology post-sputtering.

We used DFT calculations to show that S vacancies will not migrate on the surface under the experimental conditions due to the migration barriers exceeding 2 eV. This suggests that vacancy clusters form either as a result of sputtering or etching of the defective surface by oxygen. Reactions of O₂ molecules with S vacancy clusters are highly exothermic and have small or no reaction barriers. However, the kinetics of the reactions of O₂ molecule with these clusters is complex due to several factors. Our ab initio molecular dynamics simulations demonstrate that molecules are highly mobile on the surface and easily desorb at room temperature due to weak interaction. Therefore, cross-sections for reactions with vacancy dimers and trimers, which require overcoming a barrier, can be small. We demonstrate that the reaction barrier depends on the coordination of surface W atoms manifested in their effective charges – W atoms surrounded by three and more S vacancies have lower charges and are more chemically active. The O₂ dissociation on such sites proceeds without a barrier.

Another important factor is the change of the spin state of the O₂ molecule from triplet to singlet during the reaction. Although the O₂ incorporation into the WS₂ surface is accompanied by electron transfer from the surface, the reaction can be non-adiabatic, as discussed in previous studies.⁶⁶ This process requires a separate and more detailed study. Yet another factor which affects the kinetics is the fact that the reaction is highly exothermic and is accompanied by dissipation of about 8.6 eV. Modelling the mechanism of this dissipation is beyond the scope of this work. Apart from creating heat and electron–hole pairs, it may involve desorption of atoms and molecules, such as SO₂,¹⁷ which can propagate further the formation of larger vacancy clusters. As the dissociation barrier is affected by the effective charge of W atoms, this suggests that the mechanism in this work can also be extended to highly defective WS₂ monolayer. However, due to changes in the positions of the VBM and CBM between bulk and monolayer WS₂,⁷⁵ the relative position of the O₂ LUMO can influence the barrier of dissociation. Therefore, further exploration of this topic would be valuable.

In conclusion, this study provides new insights into the oxidation mechanisms of WS₂ in high defect density regimes. We show that the rate of oxidation is directly linked to the pressure of O₂ and likely the presence of exposed W atoms with low effective charges on the surface. These findings extend the understanding of WS₂ oxidation beyond the low-defect systems typically studied, offering new perspectives on how defect engineering and O₂ availability can be used to control oxidation. Future work should focus on exploring the effects of additional environmental factors, such as temperature and humidity, on the oxidation behaviour of defective WS₂, as well as strategies for passivating vacancy clusters to improve material stability.

Author contributions

D.K., F.B. and, T.D.M. performed the XPS experiments. D.K. analysed all XPS data under the supervision of F.B., M.B., A.R., and K.E.J.G. D.K. conducted all DFT and AIMD calculations under the supervision of A.S. D.K. wrote the original draft of the manuscript with contribution and editing from all authors. D.K., A.R., and A.S. edited to shape the final version of the manuscript.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

D. K. acknowledges the EPSRC Centre for Doctoral Training in Molecular Modelling & Materials Science (EP/L015862/1) and the A*STAR Graduate Academy for the graduate scholarship under the ARAP program. The authors also acknowledge the use of the ARCHER High Performance Computing Facility via membership to the UK's HPC Materials Chemistry Consortium, which is funded by EPSRC (Grant No. EP/X035859). F. B., T. D. M., and K. E. J. G. acknowledge the funding support from the Agency for Science, Technology and Research Grant (C230917006). K. E. J. G. acknowledges a Singapore National Research Foundation Grant (CRP21-2018-0001). The authors acknowledge fruitful conversations with Katherine Milton who gave insightful feedback on the manuscript and Mitisha Jain who provided insight into Ar⁺ sputtering of TMD surfaces.

References

1. W. Choi, N. Choudhary, G. H. Han, J. Park, D. Akinwande and Y. H. Lee, Mater. Today, 2017, 20, 116.
2. T. Chowdhury, E. C. Sadler and T. J. Kempa, Chem. Rev., 2020, 120, 12563–12591.
3. R. Wu, H. Zhang, H. Ma, B. Zhao, W. Li, Y. Chen, J. Liu, J. Liang, Q. Qin, W. Qi, L. Chen, J. Li, B. Li and X. Duan, Chem. Rev., 2024, 124, 10112–10191.
4. R. C. Longo, R. Addou, K. Santosh, J.-Y. Noh, C. M. Smyth, D. Barrera, C. Zhang, J. W. Hsu, R. M. Wallace and K. Cho, 2D Mater., 2017, 4, 025050.
5. Q. Li, Q. Zhou, L. Shi, Q. Chen and J. Wang, J. Mater. Chem. A, 2019, 7, 4291–4312.
6. S. Park, A. T. Garcia-Esparza, H. Abroshan, B. Abraham, J. Vinson, A. Gallo, D. Nordlund, J. Park, T. R. Kim and L. Vallez, et al. , Adv. Sci., 2021, 8, 2002768.
7. Z. Luo, W. Zheng, N. Luo, B. Liu, B. Zheng, X. Yang, D. Liang, J. Qu, H. Liu and Y. Chen, et al. , Nano Lett., 2022, 22, 2112–2119.
8. B. Chamlagain and S. I. Khondaker, ACS Omega, 2021, 6, 24075–24081.
9. M. D. Tran, J.-H. Kim and Y. H. Lee, Curr. Appl. Phys., 2016, 16, 1159–1174.
10. G. Mirabelli, C. McGeough, M. Schmidt, E. K. McCarthy, S. Monaghan, I. M. Povey, M. McCarthy, F. Gity, R. Nagle and G. Hughes, et al. , J. Appl. Phys., 2016, 120, 125102–125101.
11. A. Favron, E. Gaufrès, F. Fossard, A.-L. Phaneuf-L'Heureux, N. Y. Tang, P. L. Lévesque, A. Loiseau, R. Leonelli, S. Francoeur and R. Martel, Nat. Mater., 2015, 14, 826–832.
12. J. Gao, B. Li, J. Tan, P. Chow, T.-M. Lu and N. Koratkar, ACS Nano, 2016, 10, 2628–2635.
13. A. O. Slobodeniuk and M. R. Molas, Phys. Rev. B, 2023, 108, 035427.
14. H. Ryu, S. C. Hong, K. Kim, Y. Jung, Y. Lee, K. Lee, Y. Kim, H. Kim, K. Watanabe and T. Taniguchi, et al. , Nanoscale, 2024, 16, 5836–5844.
15. Y. Lee, J. D. S. Forte, A. Chaves, A. Kumar, T. T. Tran, Y. Kim, S. Roy, T. Taniguchi, K. Watanabe and A. Chernikov, et al. , Nat. Commun., 2021, 12, 7095.
16. L. M. Farigliano, P. A. Paredes-Olivera and E. M. Patrito, J. Phys. Chem. C, 2020, 124, 13177–13186.
17. L. M. Farigliano, P. A. Paredes-Olivera and E. M. Patrito, Phys. Chem. Chem. Phys., 2021, 23, 10225–10235.
18. Y. Guo, S. Zhou and J. Zhao, ChemNanoMat, 2020, 6, 838–849.
19. H. Liu, N. Han and J. Zhao, RSC Adv., 2015, 5, 17572–17581.
20. S. Barja, S. Refaely-Abramson, B. Schuler, D. Y. Qiu, A. Pulkin, S. Wickenburg, H. Ryu, M. M. Ugeda, C. Kastl and C. Chen, et al. , Nat. Commun., 2019, 10, 3382.
21. F. Bussolotti, H. Kawai, T. D. Maddumapatabandi, W. Fu, K. H. Khoo and K. E. J. Goh, ACS Nano, 2024, 18, 8706–8717.
22. M. Anthony and M. Seah, Surf. Interface Anal., 1984, 6, 95–106.
23. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103.
24. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488.
25. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
26. D. Kieczka, T. Durrant, K. Milton, K. E. J. Goh, M. Bosman and A. Shluger, Electron. Struct., 2023, 5, 024001.
27. J. VandeVondele and J. Hutter, Chem. Phys., 2007, 127, 114105.
28. M. Krack, Theor. Chem. Acc., 2005, 114, 145.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, Chem. Phys., 2010, 132, 154104.
30. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354.
31. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
32. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
33. R. J. Colton, A. M. Guzman and J. W. Rabalais, J. Appl. Phys., 1978, 49, 409–416.
34. C. Wagner, Faraday Discuss., 1975, 60, 291–300.
35. Y. Zhu, J. Lim, Z. Zhang, Y. Wang, S. Sarkar, H. Ramsden, Y. Li, H. Yan, D. Phuyal and N. Gauriot, et al. , ACS Nano, 2023, 17, 13545–13553.
36. F. Bussolotti, J. Yang, H. Kawai, C. P. Y. Wong and K. E. J. Goh, ACS Nano, 2021, 15, 2686–2697.
37. D. R. Baer, M. H. Engelhard, A. S. Lea, P. Nachimuthu, T. C. Droubay, J. Kim, B. Lee, C. Mathews, R. Opila and L. V. Saraf, et al. , J. Vac. Sci. Technol., B:Microelectron. Nanometer Struct. – Process., Meas., Phenom., 2010, 28, 1060–1072.
38. A. Lopez-Cazalilla, C. Cupak, M. Fellinger, F. Granberg, P. S. Szabo, A. Mutzke, K. Nordlund, F. Aumayr and R. González-Arrabal, Phys. Rev. Mater., 2022, 6, 075402.
39. J. Jussila, F. Granberg and K. Nordlund, Nucl. Mater. Energy, 2018, 17, 113–122.
40. P. Valerius, S. Kretschmer, B. V. Senkovskiy, S. Wu, J. Hall, A. Herman, N. Ehlen, M. Ghorbani-Asl, A. Grueneis and A. V. Krasheninnikov, et al. , 2D Mater., 2020, 7, 025005.
41. M. Ghorbani-Asl, S. Kretschmer, D. E. Spearot and A. V. Krasheninnikov, 2D Mater., 2017, 4, 025078.
42. Y. Shi, W. Wang, Q. Zhou, Q. Xia, D. Hua, Z. Huang, L. Chai, H. Wang and P. Wang, ACS Appl. Mater. Interfaces, 2024, 16, 29453–29465.
43. C. Kastl, R. J. Koch, C. T. Chen, J. Eichhorn, S. Ulstrup, A. Bostwick, C. Jozwiak, T. R. Kuykendall, N. J. Borys and F. M. Toma, et al. , ACS Nano, 2019, 13, 1284–1291.
44. M. P. Seah and W. Dench, Surf. Interface Anal., 1979, 1, 2–11.
45. F. Kerkhof, J. Moulijn and A. Heeres, J. Electron Spectrosc. Relat. Phenom., 1978, 14, 453–466.
46. W. Grünert, E. Shpiro, R. Feldhaus, K. Anders, G. Antoshin and K. Minachev, J. Catal., 1987, 107, 522–534.
47. A. Pashutski and M. Folman, Surf. Sci., 1989, 216, 395–408.
48. M. Kuznetsov, J. F. Zhuravlev and V. Gubanov, J. Electron Spectrosc. Relat. Phenom., 1992, 58, 169–176.
49. J. Park, I. Cho, H. Jeon, Y. Lee, J. Zhang, D. Lee, M. K. Cho, D. J. Preston, B. Shong and I. S. Kim, et al. , Adv. Mater., 2024, 36, 2314031.
50. G. Li, Y. Wang, J. Bi, X. Huang, Y. Mao, L. Luo and H. Hao, Nanomaterials, 2020, 10, 278.
51. W. Lu, B. Birmingham and Z. Zhang, Appl. Surf. Sci., 2020, 532, 147461.
52. J. Park, C. B. France and B. Parkinson, J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.–Process., Meas., Phenom., 2005, 23, 1532–1542.
53. N. S. N. Sengoku and K. O. K. Ogawa, Jpn. J. Appl. Phys., 1995, 34, 3363.
54. A. Inoue, T. Komori and K.-i. Shudo, J. Electron Spectrosc. Relat. Phenom., 2013, 189, 11–18.
55. S. Posysaev, O. Miroshnichenko, M. Alatalo, D. Le and T. S. Rahman, Comput. Mater. Sci., 2019, 161, 403–414.
56. L. Casillas-Trujillo, B. Osinger, R. Lindblad, D. Karlsson, A. I. Abrikosov, S. Fritze, K. von Fieandt, B. Alling, I. Hotz and U. Jansson, et al. , Mater. Chem. Front., 2021, 5, 5746–5759.
57. V. Nilsson, M. Van den Bossche, A. Hellman and H. Grönbeck, Surf. Sci., 2015, 640, 59–64.
58. P.-Y. Wang, C.-C. Yeh, M.-J. Chiu and C.-C. Chiu, Phys. Chem. Chem. Phys., 2024, 26, 5070–5080.
59. A. Sevy, R. F. Huffaker and M. D. Morse, J. Phys. Chem. A, 2017, 121, 9446–9457.
60. D. Tzeli, I. Karapetsas, D. M. Merriles, J. C. Ewigleben and M. D. Morse, J. Phys. Chem. A, 2022, 126, 1168–1181.
61. M. G. Bianchi, F. Risplendi, M. Re Fiorentin and G. Cicero, npj Comput. Mater., 2024, 10, 62.
62. S. Tongay, J. Zhou, C. Ataca, J. Liu, J. S. Kang, T. S. Matthews, L. You, J. Li, J. C. Grossman and J. Wu, Nano Lett., 2013, 13, 2831–2836.
63. F. R. Bagsican, A. Winchester, S. Ghosh, X. Zhang, L. Ma, M. Wang, H. Murakami, S. Talapatra, R. Vajtai and P. M. Ajayan, et al. , Sci. Rep., 2017, 7, 1774.
64. C. Zhou, W. Yang and H. Zhu, J. Chem. Phys., 2015, 142, 214706.
65. C. Schweitzer and R. Schmidt, Chem. Rev., 2003, 103, 1685–1758.
66. C. Carbogno, A. Groß, J. Meyer and K. Reuter, Dynamics of Gas-Surface Interactions: Atomic-level Understanding of Scattering Processes at Surfaces, Springer, 2013, pp. 389–419.
67. P. Gravil, J. White and D. Bird, Surf. Sci., 1996, 352, 248–252.
68. P. Gravil, D. Bird and J. White, Phys. Rev. Lett., 1996, 77, 3933.
69. M. Kurahashi, Prog. Surf. Sci., 2016, 91, 29–55.
70. D. Outka, J. Stöhr, W. Jark, P. Stevens, J. Solomon and R. Madix, Phys. Rev. B:Condens. Matter Mater. Phys., 1987, 35, 4119.
71. R. Pal, L.-M. Wang, Y. Pei, L.-S. Wang and X. C. Zeng, J. Am. Chem. Soc., 2012, 134, 9438–9445.
72. M. Baker, R. Gilmore, C. Lenardi and W. Gissler, Appl. Surf. Sci., 1999, 150, 255–262.
73. H. H. Andersen, Nucl. Instrum. Methods Phys. Res., Sect. B, 1988, 33, 466–473.
74. N. Kodama, T. Hasegawa, T. Tsuruoka, C. Joachim and M. Aono, Jpn. J. Appl. Phys., 2012, 51, 06FF07.
75. H.-g. Kim and H. J. Choi, Phys. Rev. B, 2021, 103, 085404.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4nr04992f

---