Spin-dependent vibronic response of a carbon radical ion in two-dimensional WS2

Atomic spin centers in 2D materials are a highly anticipated building block for quantum technologies. Here, we demonstrate the creation of an effective spin-1/2 system via the atomically controlled generation of magnetic carbon radical ions (CRIs) in synthetic two-dimensional transition metal dichalcogenides. Hydrogenated carbon impurities located at chalcogen sites introduced by chemical doping are activated with atomic precision by hydrogen depassivation using a scanning probe tip. In its anionic state, the carbon impurity is computed to have a magnetic moment of 1 μB resulting from an unpaired electron populating a spin-polarized in-gap orbital. We show that the CRI defect states couple to a small number of local vibrational modes. The vibronic coupling strength critically depends on the spin state and differs for monolayer and bilayer WS2. The carbon radical ion is a surface-bound atomic defect that can be selectively introduced, features a well-understood vibronic spectrum, and is charge state controlled.

F or decades, defects in semiconductors and insulators have been widely utilized to control the electronic, optical, and catalytic properties of solids. More recently, the electron spin associated with atomic defects in crystals has attracted enormous interest in light of their potential applications in quantum technologies [1][2][3] . The strong confinement of atomicscale quantum systems imposes a large characteristic energy range that is amenable for pushing quantum technologies beyond cryogenic environments 4 , while the limited system size promotes favorable coherence properties 4 .
However, generating identical defects with the necessary atomic precision, designing them to be tunable by external fields, and often even knowing their exact identities have remained unsolved challenges in the field. In this respect, two-dimensional (2D) materials offer new opportunities for the atomically-precise generation and control of defect spins 5,6 . There have been extensive efforts to utilize boron nitride or transition metal dichalcogenides (TMDs) as a platform for so-called quantum emitters, that is spin-carrying defects that can be optically detected. In boron nitride, quantum emission has been attributed to carbon impurities 7 among other defects [8][9][10] . In TMDs similar behavior has been reported but is mostly ascribed to mesoscopic strain profiles 11 . Atomic defects in TMDs have been generated previously by transmission electron microscopy 12 , and ion beam lithography 13,14 . However, the generated defects are often not identical and have mostly non-spin-polarized states.
Here, we introduce the carbon radical ion (CRI) in tungsten disulfide (WS 2 ) as an effective spin-1/2 system that can be created with atomic precision while keeping the surrounding atomic structure virtually unchanged. We precisely characterize the coupling of the spin-polarized local defect states generated by the CRI with its host lattice. Inherent to any solid-state atomic quantum system, the defect-specific electron-phonon interaction is often a dominant decay and decoherence channel 15 . Here, we demonstrate that the electron-phonon coupling associated with the two discrete electronic CRI defect states is limited to a only few vibrational modes and exhibits a distinct spin and layer dependence.
Carbon impurity defects at chalcogen sites (C X , X = S, Se) are created by scanning tunneling microscopy (STM)-induced hydrogen desorption from carbon-hydrogen (CH) complexes. Such CH impurities are frequently found in synthetic WS 2 and WSe 2 16 but can also be deliberately created by post-synthetic methane plasma treatment, as shown here for WS 2 16,17 . We induce the hydrogen depassivation of CH X by a voltage pulse from the STM tip, which is highly reproducible and atomically precise. For WS 2 , the Fermi level alignment with the graphene substrate results in a negatively charged carbon impurity with a radical character that we refer to as a CRI and denote as C À S . The CRI has an occupied spin-polarized defect state with a net magnetic moment of 1 μ B , as predicted by our density functional theory (DFT) calculations. In WSe 2 , we find C Se is neutral and thus has no magnetic moment. Moreover, we quantify the vibronic coupling of a single CRI to the host lattice by inelastic transport spectroscopy and using DFT calculations. We find that the CRI defect orbitals couple predominantly to just a few vibrational modes. However, the coupling strength critically depends on the constituent state of the CRI two-level system as well as on the number of TMD layers.

Results and discussion
In the following, we discuss our three primary conclusions: the hydrogen depassivation of the CH impurity and formation of the CRI, the two spin-polarized defect states associated with the CRI, and the vibronic coupling of the CRI defect states with the TMD host lattice (see Fig. 1).
Hydrogen depassivation of a CH impurity. In Fig. 2, STM and CO-tip noncontact atomic force microscopy (nc-AFM) images of CH À S defects in deliberately CH-doped monolayer WS 2 (0.6% atomic doping concentration) are presented. In the nc-AFM images, CH À S defects appear as a small protrusion at a sulfur site (Fig. 2c), in excellent agreement with the simulated AFM contrast obtained from the relaxed geometry using DFT calculations ( Supplementary Fig. 20). In STM, CH À S is imaged as a large, circular depression at positive bias resulting from upwards band bending due to the negative charge 18,19 . After scanning the tip over the defect at high applied sample biases and high tunneling current set-points (~2.5 V and~15 nA), a dramatic change in the STM and AFM contrast is observed. In AFM, the small protrusion disappears (see Fig. 2d). In STM, a three-fold symmetric, bright orbital structure on a dark background is observed at positive voltage. Based on the prior knowledge of the precursor defect by targeted doping, the nc-AFM simulations, and the characteristic electronic fingerprint of the converted defect (discussed below), we show that the conversion process is controlled desorption of the hydrogen atom from the CH complex. The hydrogen desorption by the STM tip is likely a resonant process where tunneling into an unoccupied CH S defect state weakens the C-H bond. This defect state as identified by DFT calculations exhibits a local anti-bonding character with a nodal plane between the carbon and hydrogen atom, supporting this hypothesis 16 . While the process is stochastic in nature, we find a minimum value of 2.3 V, required for the C-H dissociation. Tunneling at negative bias with a comparable magnitude does not result in H-desorption.
Hydrogen desorption by STM has been reported for hydrogenterminated silicon surfaces [20][21][22] and organic molecules 23,24 . As with these systems, removal of the hydrogen creates a dangling bond with a radical character, hence the defect becomes "depassivated". The dehydrogenation of the CH defect is very reliable (success rate of about 95%), and can be performed with single-atom precision as seen in Fig. 2e. Larger patterns can be Fig. 1 Creation and main attributes of the carbon radical ion (CRI). a Schematic illustration of the hydrogen depassivation of the CH impurity. b Electron transfer from the substrate into the dangling C S bond. c Magnetic moment of the carbon radical ion (CRI). d CRI vibronic coupling to the TMD host lattice. Color code: H (white), C (gray), S (yellow), W (blue), electron (small blue) with spin (red arrow).
written by scanning the surface at elevated bias and currents, shown in Fig. 2f. Occasionally, a hydrogen atom reattaches from a previously depassivated C S defect while scanning at high bias (see Supplementary Fig. 3). We speculate that an H atom residing close to the tip apex transfers back and passivates the defect again, as previously suggested 25 . This shows that the process is reversible. Even more rarely, the entire CH complex is removed irreversibly, creating a sulfur top vacancy.
Hydrogen depassivation of CH-doped TMDs can generate single defects in a 2D material with atomic precision. This is a sought-after capability for defect-based quantum systems not yet demonstrated for a 2D material. In analogy to the single transistor technologies based on hydrogen resist lithography in silicon 21,26 , we also envision that the dangling bonds of C X could be used as a predefined reactive docking site for other atoms or molecules. This approach will enable embedding functional atoms in a 2D manifold in a spatially controlled manner.
Electronic and magnetic properties of the CRI. The drastic change in STM contrast after dehydrogenation suggests a significant reconfiguration in the defect electronic structure. Scanning tunneling spectra across the CH À S defect in monolayer WS 2 are shown in Fig. 3a. The negative charge localized at the defect gives rise to a strong upwards band bending, explaining why it appears as a dark extended depression in STM images at positive bias voltage. At negative sample bias, multiple defect states are observed, which we attribute to hydrogenic bound and resonant states of the screened Coulomb potential, as we reported recently 18 .
After hydrogen desorption, two prominent defect states emerge deep in the WS 2 bandgap, one at positive (~0.6 V) and the other at negative (~−0.3 V) sample bias, as seen in Fig. 3b, d. These highly localized defect states are well decoupled from the dispersive bulk WS 2 band structure. Each state exhibits an oscillatory fine structure that is a signature of the vibronic coupling to the TMD lattice, which will be discussed later in detail. Spatial imaging of these defect states (Fig. 3f, g) reveals that they have threefold symmetry and nearly identical orbital shapes, strongly suggesting that the two resonances originate from a single open-shell defect state 27 . The gap between the defect states is slightly larger on bilayer (Δ = 770 meV) than on monolayer (Δ = 705 meV) WS 2 . Note that these numbers were corrected for the~11% voltage drop 28 across the TMD-Gr/SiC interface in our double-barrier tunneling junction geometry. The same upwards band bending before and after dehydrogenation indicates that the defect is still negatively charged, consistent with the persisting dark halo in STM images of C S .
In WSe 2 on the other hand, the carbon defect becomes chargeneutral after dehydrogenation (CH À Se → C 0 Se ). The charge neutrality can be deduced from the absence of band bending thus leading to a disappearance of the dark halo around the defect when imaged at positive sample biases upon dehydrogenation (see Supplementary Fig. 8), and the disappearance of associated hydrogenic states as seen in Fig. 3c. Moreover, the carbon impurity exhibits only a single, fully unoccupied defect state in the WSe 2 bandgap, and no state at negative bias is observed (Fig. 3e). This neutral charge is a consequence of the different band alignments of WS 2 and WSe 2 with the Gr/SiC substrate. In the WSe 2 /Gr/SiC heterostructure, the Fermi level lies roughly in the center of the WSe 2 bandgap, whereas it is 443 meV higher for WS 2 29 . Accordingly, the underlying substrate does not donate an electron to the dangling bond-like defect state of C Se in WSe 2 and thus remains empty. Alternatively, the charge state of a C X could be controlled by changing the graphene Fermi level electrostatically. 30,31 Prior DFT calculations of a neutral carbon impurity in WS 2 featured a single unoccupied defect state in the energy gap 17 . Addition of one electron to form a negatively charged carbon impurity in WS 2 results in a spin polarization and exchange-induced splitting of the in-gap state, as shown in Fig. 4b. The total energy gain of the spin-polarized ground state as compared to the less favorable non-magnetic configuration is 142 meV per single carbon atom (see Supplementary Fig. 16 for details). Hence, the negatively charged carbon impurity is characterized by a spin-split two-level system of which the lower level is occupied by one electron. Spatial maps of these defect states computed from DFT show they are strongly localized and almost identical to each other, in agreement with experimental observations ( Supplementary Fig. 9). Orbital projected density of states from our DFT calculations reveals that the defect states are strongly hybridized and possess C 2p, W 5d, and S 3p orbitals character. The charged carbon impurity is computed to possess a magnetic moment of 1 μ B from our DFT calculations, with a spin distribution closely localized at the carbon atom as shown in Fig. 4d. Based on these calculations, we conclude that the defect can be described as a CRI C À S , or CRI; and we assign the two experimentally observed in-gap states as the occupied and unoccupied spin-split defect states associated primarily with the carbon anion.
While the two defect states closely resemble a spin-1/2 system, the hybrid character of the states, with contributions from mainly  Tunneling spectroscopy of a carbon impurity in WS 2 and WSe 2 . a, b Constant height dI/dV measurement across CH À S in monolayer WS 2 before (a) and after (b) H dissociation. Both CH À S and C À S are negatively charged. The half-occupied dangling bond state of the carbon radial ion appears as two resonances in the band gap at positive and negative bias. c Constant height dI/dV measurement across C 0 Se in monolayer WSe 2 . In contrast to WS 2 , the carbon dangling bond in WSe 2 is unoccupied and the defect is neutral. d, e dI/dV spectrum recorded at the center of the C À S and C 0 Se defect, respectively. The energy difference Δ = 705 meV between electron de-and attachment in (d) results from a combination of Coulomb repulsion and spin-splitting. f, g Constant-height dI/dV maps of the C À S resonance at positive and negative bias. Their similar orbital shape indicates that they originate from the same half-occupied defect state. h Constant-height dI/dV map of the single C Se resonance. Calculated electronic and magnetic structure of a CRI. a Calculated band structure of the negatively charged CH S impurity using DFT-PBE in a 6 × 6 supercell with SOC. Note that in free-standing WS 2 the defect would be charge neutral and the valence band half-filled. Zero energy is set to the middle of the energy gap for comparison. b Calculated band structure of C À S . C S on free-standing WS 2 features only a single unoccupied defect state in the center of the WS 2 bandgap. Upon charging from the substrate, the now half-filled defect state is stabilized by spin-splitting (Δ = 593 meV, an underestimate due to the use of DFT-PBE). The gray box marks the energy range displayed in (c). c The magnetic moment of the unpaired defect spin lifts the energy degeneracy of the K and K' valley (λ = 16 meV at a defect density of 4.7 × 10 13 cm −2 ). Red and blue colors in (b) and (c) indicate opposite spin polarization. Blue (red) dots in (c) represent contribution of W l = 2, m j = −5/2(+5/2) state. d Side view of the defect-localized magnetic moment of m = 1 μ B in the out-of-plane spin configuration (for moments > 0.01 μ B ). e, f Isosurface (±0.0008 e/Å) of the difference in electron density. The occupied defect state has a higher charge density in the vicinity of the carbon atom (red). The unoccupied defect state has a higher charge density in the W plane (green). four atoms, in conjunction with SOC, leads to a small but finite magnetocrystalline anisotropy, as reported before in other solidstate spin systems 32 . Our DFT calculations predict magnetocrystalline anisotropy energy (MAE) of 0.1 meV with an easy-axis perpendicular to the WS 2 plane. Because the crystal field breaks the rotational symmetry of the impurity, the total angular momentum J is strictly speaking not a good quantum number. Rather it is more appropriate to classify defect states by the irreducible representation of the symmetry group associated with the defect site, in this case, the C 3v point group. Our DFT calculations indicate that the spatial part of both the occupied and unoccupied defect states transform symmetrically under the C 3v point group and can be labeled by the A 1 irreducible representation. Since the dimension of A 1 is one, the two defect states are non-degenerate. The Mo S antisite defect in MoS 2 with the same symmetry has been recently predicted to feature a giant MAE of about 500 meV 33 . The comparably small MAE of C À S might be expected given the small spin-orbit coupling (SOC) of carbon, similar to other light element color centers.
At ultra-low temperatures (where MAE ≫ k B T), the defect magnetic moment will be fixed in space. At a finite defect density, a net magnetic order could induce a spin-dependent shift of the spin-polarized WS 2 valence band electrons as shown in Fig. 4c, lifting the energy degeneracy between the K and K' valleys due to time-reversal symmetry breaking. This energy shift λ is calculated to be 16 meV for an ordered array of ferromagnetically coupled defects with a density of 4.7 × 10 13 cm −2 (5 × 5 unit cell) with the PBE functional 34 (see Supplementary Fig. 19 for density dependence).
The unpaired electron of the negatively charged CRI constitutes an effective spin-1/2 system. In organic chemistry, CRIs have been studied for decades in the context of reaction intermediates. Unpaired spins of organic compounds can be detected by electron spin resonance (ESR) 35 . Owing to their high reactivity, free radicals are usually very short-lived. In our experiments, the UHV conditions stabilize the CRIs, but alternatively, an inert capping layer such as hBN could be employed to protect the carbon dangling bond. It is also worth noting that isotopically pure CRIs can be easily prepared by using commercially available 12 C or 13 C clean variants of methane in the plasma treatment. Moreover, the low abundance of non-zero nuclear spin isotopes in certain TMDs and the intrinsically reduced spin densities in low-dimensional materials makes the TMD matrix a great host for defect spins 36 . While spin-bath fluctuations can be expected to act favorably in 2D TMDs, electron-phonon coupling could pose another significant source of spin decoherence, which will be discussed next.
Vibronic coupling of the CRI. Each C S defect resonance above and below the Fermi energy features characteristic, equidistant peaks, a consequence of a strong electron-phonon interaction which can be probed by the transient attachment of an electron (at positive bias voltage) or a hole (at negative bias voltage) associated with the tunneling process. Understanding the vibronic coupling of solid-state atomic quantum systems is critical as it can limit the attainable coherence times 15 . However, phonon sidebands can be effectively suppressed by the frequency-selective emission enhancement of a resonant cavity 15,37 . Moreover, lowloss local vibrations or surface acoustic phonons are analogous to nanomechanical resonators that can be used to store or coherently transmit quantum information between remote quantum systems [38][39][40] . Electron-phonon coupling is particularly relevant for very localized states, which in general lead to larger lattice relaxations. This applies to any deep center in wide-bandgap semiconductors. The theoretical framework of electron-phonon interactions at bulk defects was developed by Huang and Rhys, as well as Pekar in the 1950s and 1960s 41 . More recently, ab initio calculations of vibronic coupling to defect states in diamond 42 and in 2D materials have been a subject of active research 43,44 .
In scanning tunneling spectroscopy (STS) measurements, inelastic scattering between localized charged excitations and vibrational modes is well known to lead to characteristic sideband structures. Such phenomena have been observed for molecules on surfaces [45][46][47] , color centers in dielectrics 48 , and semiconductor quantum well tunneling devices 49 , where particular vibrational modes have been found to couple to localized electronic states. As mentioned in the previous section, when a CRI is introduced, we are able to clearly resolve this sideband structure in our dI/dV measurements (Fig. 3d). Interestingly this sideband structure differs substantially depending on the spin-state of the CRI, with the electron attachment (positive bias) exhibiting a clean Franck-Condon-like vibronic profile while the hole attachment (negative bias) exhibits a more complex fine structure, possibly involving multiple phonons. That the vibronic structure of the spin-split defect state is sensitive to the spin state is somewhat unintuitive given that these states derive from the same non-spinsplit parent state and have essentially the same orbital structure. We note that the interaction between the spin-split localized defect state and bath of harmonic phonons can be described by an effective independent spin-boson Hamiltonian 50,51 , as detailed in the Methods. From the exact solution of this model Hamiltonian, we can derive the electron spectral function A σ (ω), that is where ω ν is the frequency of the vibrational mode ν, and S νσ ¼ ðg νσ =ω ν Þ 2 is the Huang-Rhys factor, related in turn to the defect state-phonon coupling strength g νσ ; and l ν is an integer. Note that the Huang-Rhys factors and defect-phonon coupling strengths have a spin index, σ, while the vibrational frequencies do not, reflecting the fact that the frequencies are insensitive to the spin states while the coupling, in general, is not. Finally, ϵ σ denotes the electronic defect state energy, including the vibrational selfenergy. We calculate the vibrational frequencies, ω ν , using density functional perturbation theory (DFPT). Subsequently, the spindependent Huang-Rhys factors, S νσ , are extracted from a spinpolarized finite difference DFT calculation of the electronphonon coupling, using eigendisplacements from DFPT. In this manner, all quantities appearing in Eq. (5) are determined from the first principles with no adjustable parameters (for details see Methods).
In our DFT calculations of a negatively charged carbon impurity, we find two vibrational modes spatially localized around the defect that exhibit significant coupling to the defect state, with frequencies ℏω = 22 and 75 meV. The 22 meV mode is located 0.3 meV below the top of the acoustic branch of the pristine WS 2 monolayer and corresponds to an out-of-phase breathing motion involving the C-S bond and the neighboring three W atoms (Fig. 5d). Interestingly, the defect state-phonon coupling associated with this mode is very sensitive to spin and occupation of the CRI state: while we compute S = 4.5 for the unoccupied CRI state, we obtain S = 0.5 for the occupied state. This nearly order of magnitude difference can be attributed to the spin-dependence of the change in XC potential in response to specific phonon perturbations and distinct spatial parts of the wavefunctions of the two spin-polarized defect states.
To illustrate the origin of this spin-dependent coupling, we partition the Kohn-Sham Hamiltonian into two terms, a kinetic + ionic + Hartree term (T + V ext + V H ) and an XC term (V σ xc ), and express the total defect-state electron-phonon coupling as a sum of the two contributions, namely where ψ σ labels the defect state of interest and ∂ ν denotes a derivative with respect to the amplitude of displacements of phonon mode ν. The utility of the above decomposition is that it isolates the spin-dependent part of the KS Hamiltonian V σ xc . In Supplementary Fig. 13, we explicitly show how the kinetic + ionic + Hartree (KIH) term and exchange-correlation (XC) term contribute to the overall defect state-phonon coupling. For the 75 meV mode, we find for the occupied state the KIH and XC contributions add constructively while for the unoccupied state they approximately cancel out, leading to significant differences in vibrational coupling for the two states. In contrast, for the 22 meV mode, contributions from V σ xc dominate the coupling strength for the unoccupied state while it is marginal for the occupied state.
Separately, in Supplementary Fig. 13, we show the KIH term also exhibits some spin-dependence, associated with the different character of the spin-split defect state wavefunction ψ σ . To further explore this point, we plot in Fig. 4e, f the difference of the defect state electron densities, revealing the difference of the two defect states. The occupied defect state features more charge density in the vicinity of the carbon atom, which increases the coupling to the 75 meV mode, corresponding to a local out-of-plane C vibration (Fig. 5e). The unoccupied defect state has more charge density in the W plane, which leads to stronger coupling to the 22 meV mode (see Fig. 4e, f). We find that the vibronic coupling to a defect state becomes significant as the degree of the localization of the electronic wave function is sizable at the lattice sites where the vibration occurs.
From this analysis, we conclude the difference in vibrational coupling strength for the two defect states results from a combination of spin-dependent exchange and the distinct defect wavefunction character of the two states. These results are summarized in Fig. 5. As will be discussed shortly, the difference in coupling strengths is what ultimately gives rise to the very different sideband structure shown in Fig. 6a, b. Repeating our calculations for a CRI defect in a WS 2 bilayer identifies the same modes with significant coupling strengths as for monolayer, but with generally smaller S values (see Fig. 5a, b). This indicates that the defect states change less as the local vibration is excited in the bilayer system. In the bilayer, the defect states are delocalized into both layers, as shown in Supplementary  Fig. 14, reducing the coupling to the lattice vibration. For the monolayer, our calculations also reveal strong coupling to lowfrequency resonant flexural modes (Fig. 5c) that involve the CRI defect (ℏω ≈ 5 meV), particularly for the occupied CRI state. We note that a hybrid acoustic-CRI defect vibrational mode is more sensitive to supercell size, increasing the uncertainty of our calculated S values as detailed in the Supplementary Materials.
While tempting to use the DFT values for ω ν and S νσ in conjunction with Eq. (5) and compare directly with experimental STS data, we find that in practice A σ (ω) is quite sensitive to small numerical uncertainties in the Huang-Rhys factors introduced by our approximate calculations, e.g., choice of XC functional. Instead, guided by the small number of phonon modes with significant coupling strength identified by theory, we fit the measured STS spectra using Eq. (5) and subsequently compare the fit parameters to those obtained from DFT. For the unoccupied CRI state (electron attachment), three modes were sufficient for a good fit to the data; while for the occupied CRI state (hole attachment), four modes were used to largely reproduce all vibronic peaks. The fits to the tunneling spectra are shown in Fig. 6, while the refined values for ω ν and S νσ used in these fits are reported in Table 1. All vibronic resonances exhibit a homogeneous Gaussian line broadening (σ = 3-4 meV). Hence, this broadening is likely not temperature or lifetime limited but may be induced by coupling of the local vibrations to lattice acoustic modes 52 .
We pause to note that the values reported in Table 1 are not entirely unique in that different sets of ω ν and S νσ may be able to reproduce the STS data equally well. Thus care must be taken in interpreting the results. Nevertheless, there are some robust features independent of how the fit is performed which we can confidently compare with the DFT frequencies and couplings. Fig. 6 Vibronic excitations associated with charge state transitions of CRI. a Electron attachment to the unoccupied defect state for C S on monolayer WS 2 . dI/dV measurement (black dots) and three-mode electron-phonon coupling model (orange line). b Hole attachment to the occupied defect state for C S on monolayer WS 2 . dI/dV measurement (black dots) and four-mode electron-phonon coupling model (blue line). c Electron attachment for C S on bilayer WS 2 . dI/dV measurement (black dots) and three-mode electron-phonon coupling model (orange line). d Hole attachment for C S on bilayer WS 2 . dI/dV measurement (black dots) and four-mode electron-phonon coupling model (blue line). Fit parameters are given in Table 1. Notably, we find that the pristine Frank-Condon-like sideband structure observed for the unoccupied defect state (Fig. 6a) derives primarily from strong defect coupling (S ≈ 5) to a single low frequency (ℏω ≈ 20 meV) vibrational mode. Conversely, the more complex sideband structure seen in the occupied defectstate STS spectra (Fig. 6b) primarily originates from two modes, a low-frequency mode (ℏω ≈ 20 meV) with moderate coupling (S ≈ 2) and a high-frequency mode (ℏω ≈ 80 meV) with a weaker coupling (S ≈ 1), the latter manifesting a beating pattern in the STS data. Remarkably, our DFT results are consistent with these results.
In short, we identify two vibrational modes, involving the CRI defect, that couple strongly to its defect states. While all mode energies are similar for both states and layer independent, their coupling strength is greatest for monolayer WS 2 and is sensitive to the CRI defect state, indicating strong spin-phonon coupling in this system.
In summary, we demonstrate the selective and atomically precise generation of CRIs (C À S ) in a TMD host crystal. This is achieved by atomic editing of the TMD surface via STM-induced hydrogen depassivation. In its anionic state, the CRI forms an effective spin-1/2 system in the bandgap with a calculated magnetic moment of 1 μ B . Complementary measurements would be of use to confirm the calculated spin character of CRI in a magnetic field, for instance using ESR-STM. Synthetically introduced CH impurities are depassivated by H desorption using an STM tip. The resulting dangling bond introduces a deep defect state in the TMD bandgap that can be populated by electrons from the graphene substrate. For WS 2 on Gr/SiC, the Fermi level alignment leads to a negative charge state of the C S impurity, resulting in an open-shell, spin-polarized in-gap state. We also demonstrate that the atomic defect couples predominantly to two vibronic modes. While the vibrational frequencies largely defect state and layer independent, we find that the electron-phonon coupling strength is stronger for monolayer WS 2 as compared to bilayer WS 2 . The different coupling strengths to the spin-polarized defect states are a manifestation of the spin-dependent vibronic coupling in CRIs.

Methods
Sample preparation. Metalorganic chemical vapor deposition synthesis of monolayer tungsten diselenide (WSe 2 ) on graphene on silicon carbide (SiC) was performed from tungsten hexacarbonyl [W(CO) 6 , 99.99%, Sigma-Aldrich] and hydrogen selenide (H 2 Se, 99.99%, Matheson) precursors in a hydrogen gas atmosphere as previously reported. 53 Monolayer islands of tungsten disulfide (WS 2 ) on graphene/SiC were grown on graphene/SiC substrates with an ambient pressure CVD approach from tungsten oxide (WO 3 ) and sulfur powder under argon gas. 54,55 An inductively coupled plasma-enhanced CVD system was used to carry out the post-growth CH doping under an Ar/H 2 mixture and with 1 sccm methane CH 4 flow. 17 Raman spectroscopy was used to characterize the sample, as discussed previously. 16,17 After ex-situ growth, samples were transferred to ultrahigh vacuum (<2 × 10 −10 mbar) and annealed at 200°C to remove adsorbates.
Scanning probe microscopy (SPM) measurements. SPM measurements were acquired with a Createc GmbH scanning probe microscope at liquid helium temperatures (T < 7 K) under an ultrahigh vacuum (p < 2 × 10 −10 mbar). The quartz crystal cantilever (qPlus based) sensor 56 tip apex was prepared by indentations into a gold substrate and verified as metallic on the Au(111) surface. Noncontact AFM images were taken with a carbon monoxide functionalized tip 57 in constant height mode at zero bias. STM topographic measurements were taken in constant current mode with the bias applied to the sample. STS measurements were recorded using a lock-in amplifier with a resonance frequency of 683 Hz and a modulation amplitude between 2 and 10 mV.
Vibronic model. Our STS measurements of the defect states feature a characteristic sideband structure as shown in Fig. 6, arising from strong electron-phonon interactions 51,58 . In our experimental setup, resonant tunneling between the tip and sample will occur when the applied bias voltage is aligned with the defect energy levels. The defect-phonon system can be described by an effective independent spin-boson Hamiltonian 50 , that is where ϵ σ is the energy level of the defect with spin state σ, c σ ðc y σ Þ the defect annihilation (creation) operator, ω ν phonon frequency of mode ν, and a ν ða y ν Þ phonon annihilation (creation) operator, and g νσ the defect-phonon coupling matrix element.
The Hamiltonian in Eq. (3) can be solved exactly via a canonical transformation to arrive at the T = 0 electronic Green's function 50 where Φ σ ðtÞ ¼ ∑ ν S νσ ð1 À e iω ν t Þ is the phonon contribution, S νσ ¼ ðg νσ =ω ν Þ 2 the coupling constant (also known as Huang-Rhys factor), and Ω σ = ∑ ν S νσ ω ν the selfenergy. We can derive the defect energy level spectral function by taking the imaginary part of the retarded Green's function, obtaining where n denotes the number of the phonon modes that couple to the defect. In the main text, we define ϵ σ ¼ ϵ σ À Ω σ since the Ω σ only gives a relative shift of the defect level, and is not responsible for the fine structure of the sidebands. To more closely model the experimental set-up, we note that Wingreen et al. 51 added the hopping terms coupling of the contact leads with the defect state to the Hamiltonian in Eq. (3) and derived an expression for the resonant tunneling transmission probability. They found that when the hopping rate between contact and defect is assumed constant, the aforementioned transmission probability is proportional to the spectral function reported above, justifying a direct comparison of spectral function (5) with the transmission probability and ultimately the experimentally measured dI/dV curves in the main text.
First-principles calculations. In this work, the spectral function (5) is constructed in a completely ab initio manner with no adjustable parameters. We calculate ϵ σ , ω ν , and S νσ from first principles, as follows. First, for a fully optimized supercell, we perform a spin-polarized density functional theory (DFT) calculation to determine the Eigen energy associated with the defect state, ϵ σ . Using the same supercell, we then calculate ω ν at the Γ-point of the Brillouin zone using density-functional perturbation theory (DFPT) 59 . After obtaining ω ν and the corresponding (mass reduced) eigenmodes ξ ν , we construct modulated structures freezing in displacement patterns associated with all phonon modes using PHONOPY 60 . The defectphonon coupling matrix elements are then calculated with the following the finite difference method where ϵ σ (R 0 ) is the defect energy level with nuclei fixed at their equilibrium configuration, R 0 , while ϵ σ (R 0 + Aξ ν ) denotes the defect energy when nuclei are displaced, by amplitude A, along the phonon eigenmode ξ ν . For each phonon mode, we then compute the change in energy of the defect levels as a function of amplitude, A, considering a few different amplitudes. We perform a linear fit of the defect energy level with amplitude, determining g νσ from the slope [Eq. (6)]. We have explicitly confirmed, for non-spin polarized calculations, that this approach gives the same result (within 2 meV) as explicit extraction of the coupling from DFPT. In this analysis, we neglect SOC. Additional details regarding the firstprinciples calculation are provided below. Fitted Huang-Rhys factor S, vibrational energy ℏω, Gaussian broadening σ, and the electronic defect state energy ϵ of the three-mode (unoccupied state) and four-mode (occupied state) electron-phonon model shown in Fig. 6. The mode with the highest coupling strength is printed in bold.
We performed first-principles DFT calculations using the QUANTUM ESPRESSO package. 61 We employed the PBE generalized gradient approximation for XC functional 34 and used scalar relativistic optimized norm-conserving Vanderbilt (ONCV) pseudopotentials for W, S, C, and H atoms from the PseudoDojo library 62 in all cases unless explicitly including SOC within nonlinear DFT; in these cases, the fully relativistic pseudopotentials of the same type were used. We used an 80 Ry plane-wave cutoff energy. In order to model a charged system, we added one additional electron to the cell, along with compensating background charge. We used the experimental lattice parameter for WS 2 , 3.15 Å, throughout this work unless otherwise stated. We generated 3 × 3, 5 × 5, 6 × 6, and 7 × 7 supercells including a single point defect with a vacuum of 16 Å. For our total energy (electron-phonon) calculations, we used Γ-centered uniform 4 × 4(3 × 3), 3 × 3(2 × 2), 2 × 2(2 × 2), and 2 × 2(1 × 1)k-grids for 3 × 3, 5 × 5, 6 × 6, and 7 × 7 supercells, respectively, with Marzari-Vanderbilt smearing with the spread of 0.001 Ry 63 . Our convergence thresholds for ionic minimization on forces and total energy are 10 −4 and 10 −6 in atomic units. We considered spin polarization, in a collinear framework, throughout our calculations.
As the phonon frequencies of 2D materials are quite sensitive with respect to lattice parameters, we used the experimental lattice parameter of 3.15 Å since the PBE optimized lattice parameter (3.19 Å) is larger than the experimental value, which results in lower phonon frequencies. Our convergence threshold for selfconsistency in our DFPT calculations is set to 1.0 × 10 −18 a.u.

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.