# Atomic-scale characterization of defects in oxygen plasma-treated graphene by scanning tunneling microscopy

Van Dong Pham, César González, Yannick J Dappe, Chengye Dong, Joshua A Robinson, Achim Trampert, Roman Engel-Herbert

► To cite this version:

Van Dong Pham, César González, Yannick J Dappe, Chengye Dong, Joshua A Robinson, et al. Atomic-scale characterization of defects in oxygen plasma-treated graphene by scanning tunneling microscopy. Carbon, 2024, 227, 10.1016/j.carbon.2024.119260. hal-04803180

HAL Id: hal-04803180

https://cnrs.hal.science/hal-04803180v1

Submitted on 25 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Atomic-scale characterization of defects in oxygen plasma-treated graphene by scanning tunneling microscopy

Van Dong Pham<sup>a, *</sup>, César González<sup>b</sup>, Yannick J. Dappe<sup>c</sup>, Chengye Dong<sup>d</sup>,
Joshua A. Robinson<sup>d, e, f, g</sup>, Achim Trampert<sup>a</sup>, Roman Engel-Herbert<sup>a</sup>

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

Keywords:

Defect

Epitaxial graphene

SiC

Scanning tunneling microscopy/spectroscopy

Density functional theory

## ABSTRACT

Defects in graphene are important nanoscale pathways for metal atoms to enter the interface between epitaxial graphene and SiC in order to form stable ultrathin metal layers with new exotic physical properties. However, the atomic-scale details of defects that mainly govern the intercalation process remain modest. In this work, we present the first atomic investigation of point defects generated by oxygen plasma treatment on epitaxial graphene grown on SiC using low-temperature scanning tunneling microscopy, corroborated by density functional theory calculations. We found a broad spectrum of point defects that varies in size, shape, and symmetry and is dominated by triangular species. Tunneling spectroscopy identified defect-induced states in the vicinity of the Fermi level that significantly perturb the graphene electronic properties at the defect site. Based on the well-defined defect symmetry, we simulated the local density of states of the triangular defects and their corresponding scanning tunneling microscopy images which further helped us to identify the exact atomic configurations of monovacancy defects. The combination of atomic-scale scanning tunneling microscopy experiments and reliable density functional theory simulations provides ultimate microscopic details and opens a new way to identify the atomic configurations of defects in oxygen plasma-treated graphene. Our work might shed light on precise control of defect engineering in graphene for metal intercalation by controlling the defect types based on a deep understanding of each configuration.

## 1. Introduction

Point defects in graphene are important pathways to facilitate metal intercalation at the interface between epitaxial graphene (EG) and SiC [1–3]. They can be deliberately introduced into EG to promote the formation of large-scale two-dimensional (2D) metal layers which was recently established in metal confinement heteroepitaxy [4]. The basic principle of this approach is based on defect engineering controlled by oxygen plasma and thermal treatment of the graphene cap. In such context, monolayer EG is first grown on a 6H-SiC (0001) substrate by silicon thermal sublimation [5] which continuously covers the entire SiC wafer [6]. The graphene layer is then exposed to oxygen plasma to

introduce a high density of vacancy defects in the honeycomb lattice in which they act as open pathways for metal atoms to migrate into the high-energy EG/SiC interface [7] and to form crystalline ultrathin films. Hence, 2D metals are stabilized and well protected from environmental exposure thanks to the graphene cap [1,4,7] because graphene is supposed to be healed at the vacancy defects during the intercalation at a high temperature (700–800 °C) [4,8]. A broad range of experimental and theoretical efforts have attempted to access the structures of different defect types on EG and investigated how they interact with metal atoms to transport them into the EG/SiC interface [4,9,10]. These works have explicitly suggested that multivacancy defects better promote metal migration into the interface than monovacancy, thanks to

<footnote><sup>a</sup> Paul-Drude-Institut für Festkörperelektronik, Leibniz-Institut Im Forschungsverbund Berlin e. V., Hausvogteiplatz 5-7, 10117, Berlin, Germany</footnote>

<footnote><sup>b</sup> Departamento de Física de Materiales, Universidad Complutense de Madrid, E-28040, Madrid, Spain</footnote>

<footnote><sup>c</sup> SPEC, CEA, CNRS, Université Paris-Saclay, Gif-sur-Yvette Cédex, 91191, France</footnote>

<footnote><sup>d</sup> Two-dimensional Crystal Consortium, Materials Research Institute, The Pennsylvania State University, University Park, PA, 16802, United States</footnote>

<footnote><sup>e</sup> Department of Materials Science and Engineering, The Pennsylvania State University, University Park, PA, 16802, United States</footnote>

<footnote><sup>f</sup> Department of Physics, The Pennsylvania State University, University Park, PA, 16802, United States</footnote>

<footnote><sup>g</sup> Department of Chemistry, Materials Research Institute, The Pennsylvania State University, University Park, PA, 16802, United States</footnote>

<footnote>* Corresponding author.
E-mail address: pham@pdi-berlin.de (V.D. Pham).</footnote>
their large sizes and low energy barrier for metal penetration [2,3,9,11, 12].

Despite this fascinating finding which was mainly achieved based on theoretical predictions, we still lack an experimental microscopic picture allowing accurate visualization and identification of atomic configurations and electronic structures of the plasma-generated defects in graphene mainly due to the limitation in providing a direct atomic-scale characterization. Here, we use scanning tunneling microscopy (STM) and spectroscopy (STS) to provide insights into the defects in EG grown on SiC generated by oxygen plasma treatment. We found a large variety of defects characterized with atomic resolution which is mainly dominated by triangular species. This feature is representative for monovacancy and foreign atoms (i.e., oxygen) substituted monovacancy. In addition, we found other more complex defects, which exhibit an inconsistency in their topography, symmetry, size and shape that remains challenging to determine their exact atomic configurations. In corroboration with density functional theory (DFT), we simulated different atomic configurations specific for triangular defects and their corresponding density of states (DOS), allowing us to determine different triangular defect types found in the experiment. Our work has established an initial stage of atomic-scale investigation of defect-assisted metal intercalation by providing nanoscale experimental proofs of plasma-generated defects in EG on SiC using the scanning-probe technique. These results are beyond the currently available data and can potentially be exploited to further control defects in graphene and thus engineer the properties of 2D metals at the interface.

## 2. Experiment

### 2.1. Sample preparation

6H-SiC substrates (II-VI Advanced Materials) were cleaned in Acetone, IPA and nanostrip in sequence. Clean SiC was first heated up to 1500 °C and etched in 10 % H₂/Ar mixture at 700 Torr for 30 min [13]. Then monolayer EG was grown via thermal decomposition of SiC in pure argon with a pressure of 700 Torr at 1750 °C for 20 min. EG was treated with oxygen plasma in Tepla M4L plasma generator under a pressure of 0.66 bar with 150 sccm O₂ and 50 sccm He and a power of 50 W for 60s to generate defects in a controlled manner.

### 2.2. STM/STS characterization

The STM/STS measurements were conducted at sample temperature of 5 K using a cryogenic STM system (Createc) operating in ultra high vacuum (UHV) at a base pressure of 10⁻¹⁰ mbar. The STM/STS were obtained using electrochemically etched tungsten tips. The bare tip was sputtered by Ne ions and heated by electron beam under UHV condition. Before the measurements, the tips were calibrated using an Ag (111) substrate (i.e., observation of Shockley-like surface states as a step-like onset at ~ -63 mV in the tunneling spectra [14]). Alternatively, the conditioning of the tip was also achieved *in situ* on the graphene/SiC surface using field emission spectra which was described in our previous work [15]; only tips with recorded linear V-shape band of graphene and work function of about 4.2 eV typical for graphene were considered to be used in this work. The STM topography images were taken in constant-current mode. Differential conductance spectra, dI/dV versus V (where I is tunneling current and V is applied voltage) were recorded using lock-in technique with a bias modulation of 5 mV (peak-to-peak) and a modulation frequency of 675 Hz.

### 2.3. Theoretical calculations

Density Functional Theory (DFT) calculations were performed in the frame of the localized orbital basis set DFT code Fireball [16]. This code used a self-consistent version of the LDA-like Harris-Foulkes functional [17,18]. In this approach, the self-consistency of the electronic structure

calculation was achieved using the occupation numbers. Optimized numerical basis sets were used for carbon, oxygen, silicon and hydrogen, with respective cutoff radii in atomic units of s = 4.5, p = 4.5 for C [19], s = 3.3, p = 3.8 for O, s = 4.8, p = 5.4 for Si, and s = 4.1 for H. In order to describe the oxygen plasma induced defects in graphene, we considered 8 × 8 unit cells of graphene on top of a five-layer SiC slab. These big unit cells were considered to minimize the defect interactions between neighboring unit cells, due to the DFT periodicity. All the unit cells were optimized until the forces went below 0.1 eV/Å. Once the equilibrium configuration was obtained, we determined the projected density of states (PDOS) of graphene including the defects, in order to compare it with STS data. This methodology has already been applied successfully to the study of SiC van der Waals heterostructures [20,21].

For these different structures, in order to visualize and compare to the experimental STM images, we simulated the corresponding STM images from an already commonly used Keldysh-Green formalism [22, 23].

## 3. Results and discussions

First, STM imaging and tunneling spectroscopy were used to investigate the morphology and the electronic structures of EG grown on SiC, revealing important surface modification before and after oxygen plasma treatment. Fig. 1(a) shows a constant-current STM image at low sample bias (0.02 V) of pristine EG grown on SiC substrate before treatment. The graphene morphology appears as a highly crystalline monolayer free of defects. To enhance the atomic visualization, we show in the inset a close-up STM image clearly identifying the honeycomb lattice. Coexisting with the atomic resolution of EG, a larger hexagonal superstructure with a periodicity of ~18 Å is observed, indicated by the yellow rhombus. This superstructure is attributed to the 6 × 6 moiré pattern belonging to the 6√3 × 6√3 reconstruction of graphene buffer layer (GBL) indicated by the cyan cell [24–27]. This GBL is produced during the initial Si sublimation and is located right underneath the monolayer EG which has the same honeycomb structure with graphene. However, about 30 % carbon atoms of this layer are covalently bonded to the SiC substrate, having no Dirac band feature [28]. We extracted the corresponding 2D fast Fourier transform (2D-FFT) pattern from the STM image shown in Fig. 1(b), which reflects both reciprocal graphene atomic structures and the 6 × 6 moiré pattern revealed by the white and yellow hexagonal rings, respectively. Note that the 6 × 6 moiré pattern is 30° rotated with respect to the grahene orientation.

Fig. 1 (c) reveals the surface modifications of oxygen plasma-treated EG by similar surface features associated by new point defects clearly identified as bright protrusions. Notably, many defects are surrounded by an electronic scattering pattern appearing as larger hexagonal structures around each defect. This new pattern is clearly reflected in the corresponding 2D-FFT pattern shown in Fig. 1(d) as an additional purple ring. This scattering pattern is related to the quantum interferences of π-like states around the defect, giving rise to a (√3 × √3)R30 superstructure [29,30].

The electronic structures of pristine EG upon plasma treatment was measured using differential conductance (dI/dV) spectra, which is essentially proportional to the local density of state (LDOS) of the surface. Before recording the dI/dV spectra of graphene, the tip was initially calibrated on an Ag (111) substrate using repeated voltage pulses between 1 V and 3 V until the Shockley surface state appeared as an onset at ~ -63 mV in the dI/dV spectra. This indicates that the tip is purely metallic and its density of states is highly reliable for probing the LDOS of graphene. Alternatively, we also optimized the tip states by taking the dI/dV spectra in the field emission regime (FER). The emission of electrons induced by a high electrostatic field helps reform the tip apex with high stability. This approach has been successfully applied in our recent work [15]. Knowing that the Dirac point of the pristine monolayer EG on SiC occurs at about -0.4 eV below the Fermi level [31, 32], which is clearly distinguished from the phonon-induced inelastic
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Morphology and electronic properties of EG before and after oxygen plasma treatment. (a) STM topography image (0.02 V, 0.3 nA) of pristine monolayer EG/SiC. Yellow and cyan rhombuses indicate the 6 × 6 and 6√3 × 6√3 cells belonging to the GBL of the SiC substrate, respectively. (b) 2D-FFT obtained from the STM image in (a) revealing two sets of spots arranged as hexagonal rings. The outer white ring arises from the graphene 1 × 1 reciprocal lattice structure, while the inner yellow ring stems from the 6 × 6 SiC periodicity. (c) STM image (0.02 V, 0.3 nA) of monolayer EG after the oxygen plasma treatment. (d) Corresponding 2D-FFT pattern extracted from the STM image in (c) shows similar features seen in (b) with an additional purple ring emerging from the interference (√3 × √3)R30 pattern present around the bright point defects. (e) Representative dI/dV spectra measured on pristine EG (red) and on oxygen plasma-treated EG (blue) far away from defects. The black arrow marks the E<sub>Dirac</sub> around -0.4 V. (f) Distribution of the E<sub>Dirac</sub> measured on different locations on the pristine (blue bars) and the plasma treated EG (red bars) samples. (A colour version of this figure can be viewed online.)</figcaption></figure>

gap region (±67 mV) around the Fermi level, we determined the posi-
tion of the Dirac point of monolayer EG for both EG samples by
measuring the energy position of a clear minimal depth in the dI/dV
spectra recorded between –0.5 eV and 0.5 eV. It is also noted that
because at higher negative and positive bias voltages outside this range,
the dI/dV spectra often show a high instability which is not due to the
instability of the STM tip, but rather emerges from the intrinsic prop-
erties of monolayer graphene itself (see FigS. 1, Supplementary Mate-
rial). The red curve in Fig. 1(e) displays a representative dI/dV spectrum
measured on pristine EG indicating the expected linear dispersion of
both conduction and valence bands near the Dirac point (E<sub>Dirac</sub>). The
gap-like feature (~134 mV) located around the Fermi level, which
appeared as two symmetric onsets about ±67 mV is due to
phonon-mediated inelastic excitations in graphene [33,34]. The Dirac
point identified as a minimum is located at ~ –0.4 eV below the Fermi
level (marked by the black arrow). The position of the Dirac point
shifted with respect to the Fermi level clearly indicates an n-type doping
of the monolayer EG, which is well-known to be induced by the
donor-like states arising from the GBL located between the monolayer
EG and the SiC substrate [31,32] and is consistent with the previous STS
measurement [35].

The representative dI/dV spectrum of the oxygen plasma-treated
graphene taken at locations far away from defects is shown by the
blue curve in Fig. 1(e). Again, the phonon-induced gap is observed, and
the E<sub>Dirac</sub> is located at a nearly similar position as compared to that of
pristine EG. This behavior suggests that after the plasma treatment with
the introduction of new point defects into the graphene lattice, the EG is
not significantly modified in terms of charge carrier density. The sta-
tistical analysis of the Dirac energy for both graphene samples is shown
in Fig. 1(f) in order to provide a relatively precise comparison of the
charge carrier density. It first reveals a large variation in the pristine EG
(blue bars). This is mainly attributed to the inhomogeneous electronic
structure of the monolayer EG induced by the strong coupling of the SiC
substrate [27]. Note that the two dI/dV spectra shown in Fig. 1(e) only
represent the general LDOS of both graphene samples. A more detailed
picture on the graphene DOS is shown in Fig. 1(f) which reflects a large

variation in the Dirac point energy measured on different sample regions
over different experimental runs. This variation can be explained by the
fact that the monolayer EG sits directly above the SiC and is mediated by
the GBL underneath that is partly bonded to the SiC, leading to the
inhomogeneous electronic coupling between the monolayer EG and SiC
substrate on the very atomic length scale. And yet, the Dirac point is
determined from single point spectra measured at the specific atomic
sites that very likely leads to variation of the Dirac energy. It is also
worth noting that over different macroscopically distant locations where
the tip was approached, the surface impurities could be randomly pre-
sent on the graphene during the sample preparation that locally change
the EG doping behaviors. On the plasma-treated EG (red bars), the Dirac
point position is rather localized around –0.4 eV. However, it does not
reflect a clear shift, signifying no obvious modification in the charge
density of EG after the plasma treatment. Recent studies have demon-
strated that the removal of carbon atoms from graphene combined with
lattice relaxation creates a positively charged vacancy [36] (and the
references therein) which acts like an electron acceptor. In contrast, if
one assumes that the oxygen atoms which are available during the
plasma treatment substituted a certain number of empty vacancies,
these defects will be negatively charged and act like electron donor
centers [37]. The presence of both empty vacancies and
oxygen-substituted vacancies is very likely to result in an electronic
neutralization, leading to an unchanged net charge of the entire EG.
Despite that the exact amount of charge induced by the either empty or
substituted defects is not estimated, however, our experimental mean
values of the Dirac position clearly indicate no significant modification
in the net charge of EG before and after the oxygen plasma treatment. A
deeper understanding on this aspect still needs further investigation.

Our in-depth analysis of the point defect structures reveals lattice
modification with a number of different types and characteristic features
allowing a categorization. We have found a defect density of 1.4 × 10¹³
cm⁻² determined from a total of 946 defects which vary in size, shape,
apparent height and we have attempted to visualize the point defects
with atomic resolution. Fig. 2(a–c) shows three categorized triangular
defect types named as defect 1, defect 2 and defect 3. They appear not
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Atomic-scale topography and LDOS of triangular defects. (a–c) STM topography images of three distinguished triangular defects found in oxygen plasma-treated EG and classified as defect 1, defect 2, and defect 3 (at the same tunneling conditions: 0.02 V, 1 nA). Triangular defects are known to be monovacancy in the graphene lattice being either empty or substituted by foreign atoms (i.e., oxygen). (d–f) Corresponding dI/dV spectra measured above the center of defect 1, defect 2, and defect 3, respectively. Black bars highlight the positions of the localized states. (A colour version of this figure can be viewed online.)</figcaption></figure>

exactly identical, but varied slightly in the number of lobes, spatial extension from the center and their apparent height (defect 1: 1.2 Å, defect 2: 2.7 Å, defect 3: 1 Å at identical tunneling conditions). Note that the apparent height of these defects, which is mainly contributed from the density of states at the energy used to image the defects, only varies slightly when measured at different bias voltages while the atomic contrast of the graphene area around the defect is strongly bias dependent, as it can be clearly seen in FigS. 2, Supplementary Material. Three types of triangular defect were repeatedly found in different sample locations to ensure that they are explicitly originated from different atomic configurations and not due to the tip artifact. The statistical distribution of each defect type probed with atomic resolution among the three triangular defects is presented in FigS. 3, Supplementary Material, showing that the defect 3 is most frequently found.

It is known that monovacancy (MV) (which corresponds to the removal of one carbon atom) and substituted monovacancy (which corresponds to a foreign atom replacing a carbon atom) in graphene lattice, adopt a threefold symmetry as previously simulated [38,39] and experimentally identified by STM. This results from the redistribution of electron density between either the empty vacancy [40] or the substituted foreign atoms and the three nearest carbon atoms on one sublattice previously observed for nitrogen [41–44], boron [45] and hydrogen [46]. Thus, we interpret that these triangular defects are either empty MV or MV substituted by oxygen atoms considering that oxygen is the only element involved during the applied plasma treatment.

The local topographic modification as seen in the STM images implies a strong alteration in the LDOS of graphene at the defect sites unveiled by the observation of defect-induced states near Fermi level. While the spectra shown in Fig. 1(e) relate to changes of the majority carrier type in EG before and after the plasma treatment, which is probed far away from any defects, the LDOS of EG is significantly modified at the defect site itself, as shown in Fig. 2(d–f). The narrow sharp peaks indicated by black bars at the bias voltages of ~0.08 eV (defect 1), ~0.31 eV (defect 2) and ~0 eV (defect 3) are attributed to the electronic states introduced by defects, referred as defect-induced state. These sharp peaks can be interpreted by the following explanations. First, by removing one carbon atom from the graphene lattice, it creates a site with carbon dangling bonds which produces localized states. Second, the in-plane substitution of the missing carbon atoms in MV by foreign impurity (i.e. oxygen) leads to the creation of doping states which depends on the number of valence electrons of the dopant element. Here,

we believe that oxygen atoms have substituted MVs that will be further discussed in the DFT calculation presented below. Note that the position of these sharp peaks varies in the energy range around the Fermi level from defect to defect which is considered to be influenced by their local environment. This local environment around each defect is influenced by the electrostatic potential which can be easily perturbed by the interference from the nearby defects or impurities. The second possibility for this shift can be explained by further considering the relative atomic positions of the defects in the graphene lattice such as different graphene sublattices. In addition, it might also be that the inhomogeneous electronic coupling induced by the underlying GBL as revealed in Fig. 1(f) leads to the energy shift of the peaks. In either case, further detailed investigation on this aspect is remained to be further addressed. The observation of the sharp resonance peaks near Fermi level indicates that these triangular defects strongly affect the charge carrier mobility in graphene, and thus enhance the chemical reactivity at the defect sites [40] (and the references therein). It also implies that the peaks probed at the defect sites hold an important indication that the defects act like highly reactive sites for metal binding in the metal intercalation process that was theoretically predicted [4,9,10].

The modification in the LDOS of graphene induced by triangular defects is envisioned in more details in Fig. 3. It is apparent that the defect-induced state is accompanied by a significant vanishing of the phonon-induced inelastic gap at the defect site. Fig. 3(a) displays a series of consecutive dI/dV spectra measured as a function of distance when approaching a triangular defect from a clean EG area in the proximity, see Fig. 3(b). The evolution of graphene linear band and the sharp localized peak is well reflected such that at the defect center, the phonon-induced gap, highlighted by shaded area in Fig. 3(a)–is significantly vanished combined with a gradual emergence of the defect-induced state (black arrows). This trend is better detailed in the dI/dV map shown in Fig. 3(c) obtained along the dashed white line across the defect. The density on the left and right areas next to the defect decreases linearly towards the Fermi level reflecting the V-shape formed by the valence and conduction bands of the clean EG area. At the defect site, the phonon-induced inelastic gap indicated by two cyan dashed lines no longer shows a well-defined gap, and the localized state is seen as small lump within the gap, pointed by the orange arrow. In addition, along the vertical direction of the map at the defect position, the DOS is higher than that of the proximity EG area, emerging as a bright orange strip as highlighted by the yellow rectangle.

The local electronic perturbation induced by the triangular defects
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Electronic perturbation at triangular defect of EG. (a) Series of dI/dV spectra as a function of distance when going from a proximity towards the defect center shown in (b); the developments of the defect-induced state (black arrows) and the phonon-induced gap (shaded area) are clearly observed. (b) STM topography image (0.1 V, 100 pA) of a triangular defect (white circled) where the dI/dV spectra in (a) were taken. Colored dots mark the positions of the corresponding dI/dV spectra displayed in (a). (c) dI/dV map taken along the white line crossing the defect in (b). At the defect center, the phonon-induced gap (indicated by two white dashed lines) is significantly vanished together with a strong emergence of the defect-induced peak (pointed by orange arrow). The yellow dashed rectangle is used to highlight the position of the defect in the DOS map. (A colour version of this figure can be viewed online.)</figcaption></figure>

was further investigated using the spatial dI/dV map acquired with constant tip height above each defect. For instance, we measured the dI/dV map of the triangular defect 1 at the bias voltage of the resonance peak shown in Fig. 2(d) at ~0.08 V, see FigS. 4(a)(a), Supplementary Material. The triangular shape is reproduced in the map and spatially extends over several lattice sites from the defect center. This three-fold symmetry density map is originated from the charge redistribution induced either by the empty vacant or by the substituted atom in a MV [47]; the bright triangle with a higher electronic density is associated with the lower intensity extended tails facing perpendicularly to each triangular edge (see the green dashed lines in FigS. 4(a), Supplementary Material). This feature is reminiscent of STM topography of nitrogen-substituted [41–43] or boron-substituted graphene [45] in a graphitic configuration. In FigSs. 4(b) and (c), we show two other dI/dV maps obtained for defect 2 and defect 3 which were instead measured at relatively low sample biases at 0.01 V and 0.02 V, respectively and not at the exact bias voltages of the peaks found in their dI/dV spectra shown Figs. 2(e) and (f). Nevertheless, the triangular shape in the DOS maps are reproduced sharing very similar features with their STM topography images. This implies that the threefold symmetry in the DOS of triangular defects is not strictly bias-dependent in the low bias range. This can be qualitatively explained in a way that at a low sample bias, the tip integrates only the electronic states near the Fermi level (that is in most cases very few states or just a constant DOS contributing to the tunneling current) leading to a similarity between its topography and DOS map.

Knowing that the dominant triangular defects are representative for either empty or substituted monovacancies which set them apart from the other complex defects, we further proposed and theoretically investigated three different atomic structures including monovacancy (MV), an oxygen substituted carbon atom (O-MV) and three oxygens substituted carbon atoms (3O-MV) of a monovacancy using first-principle DFT. The corresponding optimized, fully relaxed structures for MV, O-MV and 3O-MV are represented in Fig. 4(a), (b) and (c), respectively. It is clearly seen from the side-view perspective of the three atomic configurations that the carbon and the substituted oxygen atoms in the graphene lattice appear atomically flat without any out-of-plane distortion, indicating no significant deformation in the lattice induced by the defects. It is important to note that the apparent height measured from the STM images shown Fig. 2(a–c) is determined by the intensity of

the tunneling current contributed by either the density of states or geometry of the surface at the defect site [48]. Therefore, it is not straightforward to interpret that these heights are real geometry of the surface at the defects. Indeed, from the relaxed structures shown in Fig. 4(a–c), one can verify that the atoms at the defect site and the neighboring carbon atoms are nearly located on the same atomic plane and hence, the defect’s height measured from the experimental STM images is contributed by the high electron density. This is because in order to maintain the constant tunneling current, the tip-sample separation should be adjusted (increased) at the defect site.

The graphene density of states was calculated including the defects to be compared with the experimental STS which reveals important proofs for identifying triangular defect types. First, the total PDOS of pristine EG was also calculated and shown in FigS. 5 (Supplementary Material) that exhibits a V-shape form as expected linear electronic band. The Dirac point is seen downshifted to −0.44 eV which is due to the electronic doping of monolayer EG induced by the SiC substrate, in good agreement with previous study [49]. The PDOS spectra calculated for three defect configurations MV, O-MV and 3O-MV are presented in Fig. 4(d), (e) and (f), respectively reveal sharp peaks located near the Fermi level. The PDOS calculated for MV exhibits two sharps peaks in Fig. 4(d). In case of O-MV, one peak around the energy range between −0.3 and −0.5 eV in Fig. 4 eV. In the PDOS calculated for 3O-MV, we find a peak located almost at the Fermi level as shown in Fig. 4(f).

After structural optimization and calculation of the electronic structure, we have simulated the corresponding STM images for MV, O-MV and 3O-MV using Keldysh-Green formalism [22,23] which are presented in Fig. 4(g), (h) and (i), respectively. At the first glance, the triangular shape is reproduced in the simulated STM images as shown in Fig. 2(a–c). Note that these images were simulated at 0.02 V which is the same sample bias obtained experimentally. The STM simulated images were also obtained at different bias voltages (see FigS. 6, Supplementary Material) in which the triangular shape typical for each defect configuration are not significantly different. This trend is consistent with the experimental bias-dependent STM images shown in Fig. 2. In the simulated STM image of the MV shown in Fig. 4(g), most of the electron density (bright lobes) is located at the three carbon atoms that are in the opposite sites (considered in the hexagonal carbon ring) of the three nearest carbon atoms next to the vacant site while the density is strongly
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Theoretical calculation of structures, density of states and STM simulation for different triangular monovacancy configurations. (a–c) Side-view and top-view of fully relaxed atomic structures showing energetically most favorable configuration of monovacancy (MV), one oxygen (O-MV) and three oxygen substituted (3O-MV) monovacancy, respectively. Grey, light yellow, red balls are C, Si, O. The pink dot in (a) highlights the missing C atoms in the MV. (d–f) PDOS of EG calculated for the three vacancies configurations, respectively, revealing localized states emerging as sharp peaks. (g–i) Corresponding simulated STM images (at 0.02 V) of the three configurations in (a), (b) and (c), respectively. (A colour version of this figure can be viewed online.)</figcaption></figure>

depleted at the vacant site. The STM simulated image for O-MV also shows a very similar feature except that the density is found localized at the three carbon atoms that are located at the opposite sites of the substituted oxygen atom in the hexagonal ring. Considering the positions of the bright lobes with respect to the carbon positions in the simulated STM images for the three cases here, one should expect a larger size for MV rather than for O-MV and 3O-MV. In case of 3O-MV, the simulated STM image shown in Fig. 4(i) show that the electron density is mainly located at the six carbon atoms adjacent to the three

substituted oxygens. This creates a distinct topographic feature as compared to the other two cases.

Combining the simulations of the three configurations shown in Fig. 4 with the experimental data in Fig. 2 allows identification of the triangular defect atomic structures. First, the simulated STM image of the 3O-MV configuration in Fig. 4 (i) convincingly correlates to the defect 3 shown in Fig. 2 (c) since they have a comparable size and the three bright lobes together with the extended tails are formed such that the positions of the three substituted oxygen atoms are clearly identified.
In addition, the center of the defect in the simulated STM image appears with low intensity, in good agreement with the experimental image. The localized state near the Fermi level as probed by the dI/dV spectrum in Fig. 2(f) and revealed by the calculated PDOS in Fig. 4(f) further supports this conclusion. Next, we attribute the O-MV shown in Fig. 4(b) to the atomic configuration of defect 2. Again, this conclusion is based on the similarity between experiment and theoretical simulations in terms of topography, spatial extensions of the three extended tails and the comparable occupied state found in the DOS. Note that the discrepancy in the DOS between experiment and simulation might arise from specific local environment around each defect measured experimentally and the size of the unit cell controlling the artificial defect-defect interaction in the calculated configurations. Among the three atomic configurations, our theoretical simulation strongly supports the identification for the cases of O-MV and 3O-MV. Regarding the MV, its PDOS and simulated STM images are partly consistent with the experimentally measured data of defect 1. On the one hand, the defect 1 appears in the experimental STM images with a quite similar size with that of defect 2 and defect 3 while the DFT simulated STM image of MV at the same bias voltage shows a considerably larger size. On the other hand, the PDOS calculated for MV displayed in Fig. 4(d) shows two strong peaks at

-0.35 eV and 0 eV, respectively, while the experimental STS spectrum of the defect 1 reveals only one peak at ~ 0.08 eV. Despite that the second peak at ~0.08 eV is reproduced and in good agreement with the calculated PDOS for MV, however, the first peak is absent in the experimental dI/dV spectra. This discrepancy can possibly be explained as the following. The experimentally measured LDOS of the defects is very sensitive to the strong electronic coupling between the defects in a monolayer EG and the SiC substrate. With a specific defect configuration, we believe that certain regions of the density of states were not clearly probed in the dI/dV spectra accounted by the tunneling matrix elements (i.e. tunneling probability at a certain bias voltage). In our opinion, this probably occurs with the defect 1. Therefore, at this stage of the measurement, we tentatively attribute the defect 1 to be MV, while great care and a deeper investigation is still needed to be accurate. Nonetheless, our simulation provides highly reliable approach to identify the atomic and electronic structures of triangular defects, which establishes principle for investigating more complex structures found in our sample.

Apart from the triangular defects, we found a broad spectrum of other defects with extended size and shape, and different symmetry compared to the triangular defects. As shown in the statistical graph of

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Exemplified complex defects found in addition with the triangular defects. (a-d, left column) STM topography images of few example complex defect found in the oxygen plasma-treated EG sample in addition to the triangular monovacancy defects. (e-h, middle column) Corresponding dI/dV spectra measured above each defect in the left column; the sharp peaks are indicated by black bars. (i-l, right column) dI/dV maps acquired at the energies either at the peak positions of the dI/dV spectra shown in middle column, or at low bias voltage. (A colour version of this figure can be viewed online.)</figcaption></figure>
defect types (Fig. S 3, Supplementary Material) determined from total 92
defects imaged with atomic resolution, about 35 % of these defects are
not triangular suggesting that more complex structures with atomic
configurations beyond monovacancy (i.e. multivacancy) were intro-
duced into graphene by the oxygen plasma treatment.

In Fig. 5 (a-d, left column) we exemplify the topography of some
complex defects with high atomic resolution found in the plasma-treated
EG sample. It is apparent that these defects are not identical and do not
appear with a highly rotational symmetry. Their corresponding dI/dV
spectra shown in Fig. 5(e-h, middle column) reveal defect-induced states
as narrow peaks located in either valence or conduction bands near the
Fermi level. In the right column of Fig. 5 (i-l) we show dI/dV maps of
these defects acquired either at low bias voltages or at the precise peak
energy of the corresponding dI/dV spectra. The DOS maps mainly
replicate their topography STM images. The narrow peaks in the dI/dV
spectroscopy and the local DOS map measured on these complex defects
show quite similar tendency to those measured for the triangular de-
fects. This again implies an enhancing chemical reactivity at these sites
with adsorbates. This has been demonstrated in case of nitrogen-doped
EG on SiC, in which the nitrogen doping sites give rise to a strong charge
transfer between the dopant sites of graphene and the physisorbed
organic molecules [50,51]. In such context, it should be straightforward
to expect a strong interaction between metal atoms and the defect sites
when using defects as major pathways for metal intercalation.

Previous DFT studies have predicted multivacancies in the treated
EG taking into account Jahn-Teller distortion bond reconstruction [9]. A
broad spectrum of bare vacancy ranging from mono to octa structures
were simulated. In comparison with experimental data found here, we
do not find evidences that are directly relevant to the configurations
described in that work. This is because even though STM images do not
allow a direct identification of atom-atom bonding and structures,
nevertheless, the symmetry and the atomic contour of the multivacancy
should be partially reflected in the topography of a specific configura-
tion. The defects listed in Fig. 5 however appear highly random sug-
gesting that they adopt even more complex structures. In our best
understanding, this can be attributed to the random removal of carbon
atoms from the honeycomb lattice under plasma treatment. And in
addition to the creation of multivacancies, the asymmetric features of
the defects shown in Fig. 5 may also suggest different origin: that is the
deformation of the carbon lattice leading to wrinkle topographic fea-
tures around the defects under plasma bombardment or out-of-plane
oxygen (or other unknown impurities) decorated graphene at the
vacant sites. In any case, it still remains challenging in identifying the
atomic configurations of complex defects based only on the STM
topography and dI/dV spectroscopy, since there is no indication of
which atomic configuration should be theoretically simulated, hence
further works are still remained to be investigated. Nevertheless, the
close-up inspection at atomic level in terms of topography, electronic
properties and statistics on types of defects in our work has provided a
convincing atomic picture of the oxygen plasma-treated epitaxial gra-
phene on SiC.

4. Conclusions

In summary, we have investigated the structures and electronic
properties of different defect types on epitaxial graphene grown on SiC
which are generated by oxygen plasma treatment. Atomic resolution
STM topography and tunneling spectroscopy allow us to access the
atomic-scale details of defects ranging from triangular shape (that are
monovacancy and its related species) to more complex structures in
corroboration with first-principle DFT calculations. From the density of
state viewpoint, the defect-induced states are localized near the Fermi
level and significantly modify the band structure of EG at the defect sites
themselves. This observation indicates that the chemical reactivity of
graphene is significantly enhanced at the defect sites, reflecting their
important role in facilitating the metal migration in the EG/SiC

interface. Dominant triangular defects with well-defined symmetry
found in the STM measurement allow simulation of their topographies
and DOS from which their atomic configurations were precisely iden-
tified. Despite the atomic configurations of other complex defects have
not yet been corroborated, nonetheless detailed measurements of their
topographies and LDOS provide important experimental proofs which
advance the future investigations to be carried out. The significance of
our work is perhaps best realized by deeply analyzing the atomic
structures of defects and their electronic properties at atomic level and
may open new way for visualizing the metal intercalation mechanism in
which this picture at nanometer scale will be explored using scanning-
probe approach.

CRediT authorship contribution statement

**Van Dong Pham:** Writing – review & editing, Writing – original draft, Investigation, Conceptualization. **César González:** Investigation. **Yannick J. Dappe:** Investigation. **Chengye Dong:** Investigation, Conceptualization. **Joshua A. Robinson:** Investigation, Conceptualization. **Achim Trampert:** Writing – review & editing, Investigation. **Roman Engel-Herbert:** Investigation.

Declaration of competing interest

The authors declare the following financial interests/personal re-
lationships which may be considered as potential competing interests:

Joshua A. Robinson, Chengye Dong reports financial support was
provided by 2DCC-MIP under NSF cooperative agreement DMR-
2039351 and the Center for Nanoscale Science through the NSF Grant
DMR-2011839. If there are other authors, they declare that they have no
known competing financial interests or personal relationships that could
have appeared to influence the work reported in this paper.

Acknowledgements

The authors thank J.M.J. Lopes (Paul-Drude-Institut für Festköper-
elektronik), V.-T. Tran (Université Paris-Saclay, CNRS, France) for crit-
ical reading the manuscript and S. Fölsch (Paul-Drude-Institut für
Festköperelektronik) for technical supports in this work. J. A. Robinson
and C. Dong are supported by 2DCC-MIP under NSF cooperative
agreement DMR-2039351 and the Center for Nanoscale Science through
the NSF Grant DMR-2011839.

**Appendix A. Supplementary data**

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.carbon.2024.119260.

**References**

[1] N. Briggs, Z.M. Gebeyehu, A. Vera, T. Zhao, K. Wang, A. De La Fuente Duran, B. Bersch, T. Bowen, K.L. Knappenberger, J.A. Robinson, Epitaxial graphene/ silicon carbide intercalation: a minireview on graphene modulation and unique 2D materials, Nanoscale 11 (2019) 15440–15447, https://doi.org/10.1039/ C9NR03721G.

[2] X. Lu, Y. Liu, M. Shao, X. Liu, Defect-mediated intercalation of dysprosium on buffer layer graphene supported by SiC(0001) substrate, Chem. Phys. Lett. 742 (2020), https://doi.org/10.1016/j.cplett.2020.137162.

[3] Y. Liu, X. Liu, C.-Z. Wang, Y. Han, J.W. Evans, A. Lii-Rosales, M.C. Tringides, P. A. Thiel, Mechanism of metal intercalation under graphene through small vacancy defects, J. Phys. Chem. C 125 (2021) 6954–6962, https://doi.org/10.1021/acs. jpcx.1c00814.

[4] N. Briggs, B. Bersch, Y. Wang, J. Jiang, R.J. Koch, N. Nayir, K. Wang, M. Kolmer, W. Ko, A. De La Fuente Duran, S. Subramanian, C. Dong, J. Shallenberger, M. Fu, Q. Zou, Y.-W. Chuang, Z. Gai, A.-P. Li, A. Bostwick, C. Jozwiak, C.-Z. Chang, E. Rotenberg, J. Zhu, A.C.T. Van Duin, V. Crespi, J.A. Robinson, Atomically thin half-van der Waals metals enabled by confinement heteroepitaxy, Nat. Mater. 19 (2020) 637–643, https://doi.org/10.1038/s41563-020-0631-x.

[5] C. Berger, Z. Song, T. Li, X. Li, A.Y. Ogbazghi, R. Feng, Z. Dai, A.N. Marchenkov, E. H. Conrad, P.N. First, W.A. De Heer, Ultrathin epitaxial graphite: 2D electron gas
properties and a route toward graphene-based nanoelectronics, J. Phys. Chem. B 108 (2004) 19912–19916, https://doi.org/10.1021/jp040650f.

[6] C. Riedl, C. Coletti, T. Iwasaki, A.A. Zakharov, U. Starke, Quasi-free-standing epitaxial graphene on SiC obtained by hydrogen intercalation, Phys. Rev. Lett. 103 (2009) 246804, https://doi.org/10.1103/PhysRevLett.103.246804.

[7] C. Berger, W.A. de Heer, Flat and safe under the graphene sheet, Nat. Mater. 19 (2020) 583–584, https://doi.org/10.1021/s41563-020-0666-z.

[8] J. Chen, T. Shi, T. Cai, T. Xu, L. Sun, X. Wu, D. Yu, Self healing of defected graphene, Appl. Phys. Lett. 102 (2013) 103107, https://doi.org/10.1063/1.4795292.

[9] N. Nayir, M.Y. Sengul, A.L. Costine, P. Reinke, S. Rajabpour, A. Bansal, A. Kozhakhmetov, J. Robinson, J.M. Redwing, A. Van Duin, Atomic-scale probing of defect-assisted Ga intercalation through graphene using ReaxFF molecular dynamics simulations, Carbon 190 (2022) 276–290, https://doi.org/10.1016/j.carbon.2022.01.005.

[10] F. Niefind, Q. Mao, N. Nayir, M. Kowalik, J. Ahn, A.J. Winchester, C. Dong, R. A. Maniyara, J.A. Robinson, A.C.T. Van Duin, S. Pookpanratana, Watching (De) Intercalation of 2D metals in epitaxial graphene: insight into the role of defects, Small (2023) 2306554, https://doi.org/10.1002/smll.202306554.

[11] Y. Orimoto, K. Otsuka, K. Yagyu, H. Tochihara, T. Suzuki, Y. Aoki, Theoretical study of Cu intercalation through a defect in zero-layer graphene on SiC surface, J. Phys. Chem. C 121 (2017) 7294–7302, https://doi.org/10.1021/acs.jpcb.7b00314.

[12] Y. Zhou, A. Lii-Rosales, M. Kim, M. Wallingford, D. Jing, M.C. Tringides, C.-Z. Wang, P.A. Thiel, Defect-mediated, thermally-activated encapsulation of metals at the surface of graphite, Carbon 127 (2018) 305–311, https://doi.org/10.1016/j.carbon.2017.10.103.

[13] S. Subramanian, D.D. Deng, K. Xu, N. Simonson, K. Wang, K. Zhang, J. Li, R. Feenstra, S.K. Fullerton-Shirey, J.A. Robinson, Properties of synthetic epitaxial graphene/molybdenum disulfide lateral heterostructures, Carbon 125 (2017) 551–556, https://doi.org/10.1016/j.carbon.2017.09.058.

[14] M. Pivetta, F. Silly, F. Patthey, J.P. Pelz, W.-D. Schneider, Reading the ripples of confined surface-state electrons: profiles of constant integrated local density of states, Phys. Rev. B 67 (2003) 193402, https://doi.org/10.1103/PhysRevB.67.193402.

[15] V.D. Pham, C. Dong, J.A. Robinson, Atomic structures and interfacial engineering of ultrathin indium intercalated between graphene and a SiC substrate, Nanoscale Adv. 5 (2023) 5601–5612, https://doi.org/10.1039/D3NA00630A.

[16] J.P. Lewis, P. Jelínek, J. Ortega, A.A. Demkov, D.G. Trabada, B. Haycock, H. Wang, G. Adams, J.K. Tomfohr, E. Abad, H. Wang, D.A. Drabold, Advances and applications in the FIREBALL ab initio tight-binding molecular-dynamics formalism, Phys. Status Solidi B 248 (2011) 1989–2007, https://doi.org/10.1002/pssb.201147259.

[17] J. Harris, Simplified method for calculating the energy of weakly interacting fragments, Phys. Rev. B 31 (1985) 1770–1779, https://doi.org/10.1103/PhysRevB.31.1770.

[18] W.M.C. Foulkes, R. Haydock, Tight-binding models and density-functional theory, Phys. Rev. B 39 (1989) 12520–12536, https://doi.org/10.1103/PhysRevB.39.12520.

[19] M.A. Basanta, Y.J. Dappe, P. Jelínek, J. Ortega, Optimized atomic-like orbitals for first-principles tight-binding molecular dynamics, Comput. Mater. Sci. 39 (2007) 759–766, https://doi.org/10.1016/j.commatsci.2006.09.003.

[20] M.T. Dau, M. Gay, D. Di Felice, C. Vergnaud, A. Marty, C. Beigné, G. Renaud, O. Renault, P. Mallet, T. Le Quang, J.-Y. Veuilllen, L. Huder, V.T. Renard, C. Chapelier, G. Zamborlini, M. Jugovac, V. Feyer, Y.J. Dappe, P. Pochet, M. Jamet, Beyond of the Waals Interaction: The Case of MoSe₂ Epitaxially Grown on Few-Layer Graphene, ACS Nano 12 (2018) 2319–2331, https://doi.org/10.1021/acsnano.7b07446.

[21] Y.J. Dappe, Y. Almadori, M.T. Dau, C. Vergnaud, M. Jamet, C. Paillet, T. Journot, B. Hyot, P. Pochet, B. Grévin, Charge transfers and charged defects in WSe₂ graphene-SiC interfaces, Nanotechnology 31 (2020) 255709, https://doi.org/10.1088/1361-6528/ab8083.

[22] C. González, E. Abad, Y.J. Dappe, J.C. Cuevas, Theoretical study of carbon-based tips for scanning tunneling microscopy, Nanotechnology 27 (2016) 105201, https://doi.org/10.1088/0957-4484/27/10/105201.

[23] S. Dubey, S. Lisi, G. Nayak, F. Herziger, V.-D. Nguyen, T. Le Quang, V. Cherkez, C. González, Y.J. Dappe, K. Watanabe, T. Taniguchi, L. Magaud, P. Mallet, J.-Y. Veuillen, R. Arenal, L. Marty, J. Renard, N. Bendiab, J. Coraux, V. Bouchiat, Weakly trapped, charged, and free excitons in single-layer MoS₂ in the presence of defects, strain, and charged impurities, ACS Nano 11 (2017) 11206–11216, https://doi.org/10.1021/acsnano.7b05520.

[24] G.M. Rutter, N.P. Guisinger, J.N. Crain, E.A.A. Jarvis, M.D. Stiles, T. Li, P.N. First, J.A. Stroscio, Imaging the interface of epitaxial graphene with silicon carbide via scanning tunneling microscopy, Phys. Rev. B 76 (2007) 235416, https://doi.org/10.1103/PhysRevB.76.235416.

[25] C. Riedl, C. Coletti, U. Starke, Structural and electronic properties of epitaxial graphene on SiC(0 0 0 1): a review of growth, characterization, transfer doping and hydrogen intercalation, J. Phys. D Appl. Phys. 43 (2010) 374009, https://doi.org/10.1088/0022-3727/43/37/374009.

[26] T.W. Hu, D.Y. Ma, F. Ma, K.W. Xu, Preferred armchair edges of epitaxial graphene on 6H-SiC(0001) by thermal decomposition, Appl. Phys. Lett. 101 (2012) 241903, https://doi.org/10.1063/1.4769967.

[27] V.W. Brar, Y. Zhang, Y. Yayon, T. Ohta, J.L. McChesney, A. Bostwick, E. Rotenberg, K. Horn, M.F. Crommie, Scanning tunneling spectroscopy of inhomogeneous electronic structure in monolayer and bilayer graphene on SiC, Appl. Phys. Lett. 91 (2007) 122102, https://doi.org/10.1063/1.2771084.

[28] F. Varchon, R. Feng, J. Hass, X. Li, B.N. Nguyen, C. Naud, P. Mallet, J.-Y. Veuillen, C. Berger, E.H. Conrad, L. Magaud, Electronic structure of epitaxial graphene layers on SiC: effect of the substrate, Phys. Rev. Lett. 99 (2007) 126805, https://doi.org/10.1103/PhysRevLett.99.126805.

[29] G.M. Rutter, J.N. Crain, N.P. Guisinger, T. Li, P.N. First, J.A. Stroscio, Scattering and interference in epitaxial graphene, Science 317 (2007) 219–222, https://doi.org/10.1126/science.1142882.

[30] Y. Zhang, Y. Su, L. He, Local berry phase signatures of bilayer graphene in intervalley quantum interference, Phys. Rev. Lett. 125 (2020) 116804, https://doi.org/10.1103/PhysRevLett.125.116804.

[31] Y. Qi, S.H. Rhim, G.F. Sun, M. Weinert, L. Li, Epitaxial graphene on SiC(0001): more than just honeycombs, Phys. Rev. Lett. 105 (2010) 085502, https://doi.org/10.1103/PhysRevLett.105.085502.

[32] J. Ristein, S. Mammadov, Th Seyller, Origin of doping in quasi-free-standing graphene on silicon carbide, Phys. Rev. Lett. 108 (2012) 246104, https://doi.org/10.1103/PhysRevLett.108.246104.

[33] Y. Zhang, V.W. Brar, F. Wang, C. Girit, Y. Yayon, M. Panlasigui, A. Zettl, M. F. Crommie, Giant phonon-induced conductance in scanning tunnelling spectroscopy of gate-tunable graphene, Nat. Phys. 4 (2008) 627–630, https://doi.org/10.1038/nphys1022.

[34] V.W. Brar, S. Wickenburg, M. Panlasigui, C.-H. Park, T.O. Wehling, Y. Zhang, R. Decker, Ç. Girit, A.V. Balatsky, S.G. Louie, A. Zettl, M.F. Crommie, Observation of carrier-density-dependent many-body effects in graphene via tunneling spectroscopy, Phys. Rev. Lett. 104 (2010) 036805, https://doi.org/10.1103/PhysRevLett.104.036805.

[35] Y. Liu, M. Weinert, L. Li, Determining charge state of graphene vacancy by noncontact atomic force microscopy and first-principles calculations, Nanotechnology 26 (2015) 035702, https://doi.org/10.1088/0957-4484/26/3/035702.

[36] J. Mao, Y. Jiang, D. Moldovan, G. Li, K. Watanabe, T. Taniguchi, M.R. Masir, F. M. Peeters, E.Y. Andrei, Realization of a tunable artificial atom at a supercritically charged vacancy in graphene, Nat. Phys. 12 (2016) 545–549, https://doi.org/10.1038/nphys3665.

[37] D.M.A. Mackenzie, M. Galbiati, X.D. De Cerio, I.Y. Sahalianov, T.M. Radchenko, J. Sun, D. Peña, L. Gammelgaard, B.S. Jessen, J.D. Thomsen, P. Bøggild, A. Garcia-Lekue, L. Camilli, J.M. Caridad, Unraveling the electronic properties of graphene with substitutional oxygen, 2D Mater. 8 (2021) 045035, https://doi.org/10.1088/2053-1583/ac28ab.

[38] H.A. Mizes, J.S. Foster, Long-range electronic perturbations caused by defects using scanning tunneling microscopy, Science 244 (1989) 559–562, https://doi.org/10.1126/science.244.4904.559.

[39] Ph Lambin, H. Amara, F. Ducastelle, L. Henrad, Long-range interactions between substitutional nitrogen dopants in graphene: electronic properties calculations, Phys. Rev. B 86 (2012) 045448, https://doi.org/10.1103/PhysRevB.86.045448.

[40] M.M. Ugeda, I. Brihuega, F. Guinea, J.M. Gómez-Rodríguez, Missing atom as a source of carbon magnetism, Phys. Rev. Lett. 104 (2010) 096804, https://doi.org/10.1103/PhysRevLett.104.096804.

[41] L. Zhao, R. He, K.T. Rim, T. Schiros, K.S. Kim, H. Zhou, C. Gutiérrez, S. P. Chockalingam, C.J. Arguello, L. Pálová, D. Nordlund, M.S. Hybertsen, D. R. Reichman, T.F. Heinz, P. Kim, A. Pinczuk, G.W. Flynn, A.N. Pasupathy, Visualizing individual nitrogen dopants in monolayer graphene, Science 333 (2011) 999–1003, https://doi.org/10.1126/science.1208759.

[42] F. Joucken, Y. Tison, J. Lagoute, J. Dumont, D. Cabosart, B. Zheng, V. Repain, C. Chacon, Y. Girard, A.R. Botello-Méndez, S. Rousset, F. Sporken, J.-C. Charlier, L. Henrad, Localized state and charge transfer in nitrogen-doped graphene, Phys. Rev. B 85 (2012) 161408, https://doi.org/10.1103/PhysRevB.85.161408.

[43] F. Joucken, L. Henrad, J. Lagoute, Electronic properties of chemically doped graphene, Phys. Rev. Mater. 3 (2019) 110301, https://doi.org/10.1103/PhysRevMaterials.3.110301.

[44] J. Lagoute, F. Joucken, V. Repain, Y. Tison, C. Chacon, A. Bellec, Y. Girard, R. Sporken, E.H. Conrad, F. Ducastelle, M. Palsgaard, N.P. Andersen, M. Brandbyge, S. Rousset, Giant tunnel-electron injection in nitrogen-doped graphene, Phys. Rev. B 91 (2015) 125442, https://doi.org/10.1103/PhysRevB.91.125442.

[45] L. Zhao, M. Levendorf, S. Goncher, T. Schiros, L. Pálová, A. Zabet-Khosousi, K. T. Rim, C. Gutiérrez, D. Nordlund, C. Jaye, M. Hybertsen, D. Reichman, G. W. Flynn, J. Park, A.N. Pasupathy, Local atomic and electronic structure of boron chemical doping in monolayer graphene, Nano Lett. 13 (2013) 4659–4665, https://doi.org/10.1021/nl401781d.

[46] D. Wong, Y. Wang, W. Jin, H.-Z. Tsai, A. Bostwick, E. Rotenberg, R.K. Kawakami, A. Zettl, A.A. Mostofi, J. Lischner, M.F. Crommie, Microscopy of hydrogen and hydrogen-vacancy defect structures on graphene devices, Phys. Rev. B 98 (2018) 155436, https://doi.org/10.1103/PhysRevB.98.155436.

[47] K.F. Kelly, D. Sarkar, G.D. Hale, S.J. Oldenburg, N.J. Halas, Threefold electron scattering on graphite observed with C 60 -adsorbed STM tips, Science 273 (1996) 1371–1373, https://doi.org/10.1126/science.273.5280.1371.

[48] C. Park, M. Yoon, Topography inversion in scanning tunneling microscopy of single-atom-thick materials from penetrating substrate states, Sci. Rep. 12 (2022) 7321, https://doi.org/10.1038/s41598-022-10870-0.
[49] A. Bostwick, T. Ohta, T. Seyller, K. Horn, E. Rotenberg, Quasiparticle dynamics in graphene, Nat. Phys. 3 (2007) 36–40, https://doi.org/10.1038/nphys477.

[50] V.D. Pham, J. Lagoute, O. Mouhoub, F. Joucken, V. Repain, C. Chacon, A. Bellec, Y. Girard, S. Rousset, Electronic interaction between nitrogen-doped graphene and porphyrin molecules, ACS Nano 8 (2014) 9403–9409, https://doi.org/10.1021/nn503753e.

[51] V.D. Pham, F. Joucken, V. Repain, C. Chacon, A. Bellec, Y. Girard, S. Rousset, R. Sporken, M.C.D. Santos, J. Lagoute, Molecular adsorbates as probes of the local properties of doped graphene, Sci. Rep. 6 (2016) 24796, https://doi.org/10.1038/srep24796.