Sample preparation<\/p>\n
Our WSe2 samples are grown via metal organic chemical vapor deposition on few-layer epitaxial graphene supported by a conductive silicon carbide (c-SiC) substrate. The plasma frequency of c-SiC is sufficient to effectively suppress THz reflections from the bottom surface. Intercalation of the graphene layers with hydrogen prior to WSe2 growth terminates the dangling silicon bonds, providing a homogeneous electrostatic background. Further details of the sample growth procedure were previously reported along references16<\/a>,43<\/a>. Key advantages of epitaxially grown samples on quasi-free-standing graphene (QFEG) are the electronic homogeneity (no moir\u00e9), a Fermi energy pinned in the center of the large WSe2 band gap, and weak hybridization with the substrate16<\/a>. The strong out-of-plane confinement facilitates electron-phonon coupling, well-defined defect states, and high contrast of the defect states inside the band gap, enabling LW-STM at nonzero Vdc.<\/p>\n High-temperature annealing in UHV is required to clean samples after ambient exposure. We generate isolated VacSe in the top WSe2 layer via light sputtering with Ar ions, approximately 2\u2009s at 0.12\u2009kV acceleration bias, 0.16\u2009kV discharge voltage, and a sub-\u03bcA sputter-current16<\/a>. Multiple VacSe in 1\u20132\u2009ML WSe2 were investigated for this study, emphasizing the robustness and repeatability of the technique. Fig.\u2009S1<\/a> lists STM topography, dI\/dV spectra, and current saturation curves of the relevant VacSe.<\/p>\n STM height set point and average charge-state lifetime<\/p>\n For a reproducible definition of the tip height z0, the STM feedback loop was opened on pristine WSe2 at a set point current of 100\u2009pA with applied bias voltage of 1.5\u2009V (1\u2009ML) or 1.2\u2009V (2\u2009ML), corresponding to an energy 200\u2009meV above the conduction band onset. We estimate an absolute tip\u2013sample distance at z0 of (7\u2009\u00b1\u20092)\u2009\u00c5. The average charge-state lifetime (\u03c40) of VacSe states scales with the decoupling layers of the sample and is determined via the current saturation in STM approach curves15<\/a>,16<\/a>,32<\/a>. For simplicity, we consider only the first unoccupied defect state (LUMO) of VacSe for the determination of \u03c40. Reduced localization due to a lower binding energy slightly reduces the average charge-state lifetime of the LUMO+1 defect state16<\/a>. For different VacSe, \u03c40 slightly vary due to differences in their dielectric environment, however, the energy spectra and orbital shape remain comparable within a few-10\u2009meV shifts (Fig.\u2009S1<\/a>).<\/p>\n Voltage drop and definition of \u0394V<\/p>\n In our double barrier tunneling junction geometry, the effective voltage from the tip to the VacSe defect states depends on the WSe2 dielectric screening and the absolute height set point z0. The relative voltage \u0394V\u2009=\u2009Vdc\u2009\u2212\u2009VLUMO(z) is important to compare LW-STM data across different VacSe, where the onset of the LUMO, VLUMO(z), varies on the order of tens of meV due to inhomogeneity of the sample and as a function of z. While for \u0394V\u2009<\u20090 the dc current is negligible, \u0394V\u2009>\u20090 V generates a notable dc current and thus increases the\u00a0VacSe population\u00a0statically in the experiment. In orbital imaging, we define \u0394V at a single point on the orbital lobe of VacSe, avoiding complications of lateral dependence. We determine VLUMO(z) directly from a dI\/dV measurement at z or via extrapolation of VLUMO(z0) on the basis of measurements performed in ref. 16<\/a>, Fig.\u20092<\/a>. Owing to the high sensitivity towards thermal drifts on the sub-\u00c5 scale, we estimate an accuracy \u03b4\u0394V\u2009\u2248\u200910\u2009mV.<\/p>\n LW-STM setup and near-field waveform detection<\/p>\n Our LW-STM system is a commercial low-temperature STM from CreaTec Fischer & Co. GmbH with free-space optical access to the STM junction via large numerical aperture parabolic mirrors focused on the tip apex. THz pulses are generated via tilted pulse front optical rectification in lithium niobate pumped by a multi-MHz Yb:fiber laser with pump pulse energies of a few \u03bcJ. The system achieves peak THz voltages at the STM tip up to 0.5, 1, and 2\u2009V at 41, 20, and 10\u2009MHz repetition rates, respectively. We use frustrated internal reflection in a PTFE prism as well as a double or single mirror sequence to adapt the phase and polarity of single-cycle THz pulses at the STM tip44<\/a>. Pump and probe pulses are generated independently and are\u00a0recombined in collinear geometry via a c-cut sapphire window. Purging of the optics setup with dry air reduces ambient absorption before pulses are coupled into the vacuum chamber. Operating conditions of the scanning probe microscope are 5\u2009K base temperature and 10\u221210\u2009mbar pressure.<\/p>\n To discriminate lightwave-driven currents, we mechanically modulate the THz path at an intermediate focus using frequencies between 800 and 1000\u2009Hz that are compatible with the 109 transimpedance amplifier of the STM and show weak mechanical and electric noise. ILW is then extracted via lock-in demodulation and calibrated with respect to the absolute current. We convert current to rectified charge per transient (e\/tr) on the basis of the laser repetition rate and the 50% duty cycle of the modulation. In pump-probe measurements, only the probe pulse is modulated, while pump pulses enter the junction at full repetition rate. We find that high DC currents prevent sensitive measurements of QLW because of large noise floors.<\/p>\n The THz peak voltage is calibrated against Vdc using a recently developed approach based on the strongly nonlinear onset of the WSe2 conduction band31<\/a> measured in a clean pristine region of the sample (Fig.\u2009S2a\u2013c<\/a>). On the basis of this amplitude calibration, we determine the relative time origin (\u0394t\u2009=\u20090) and temporal resolution of the experiment (380\u2009fs) by cross correlation of equal THz fields where only the sum field rectifies charges (Fig.\u2009S2d<\/a>). In a second step, we measure the transient near-field waveform at the STM tip using THz cross correlation (THz-CC)29<\/a>,31<\/a>. Here, we set relative amplitudes such that a strong THz gate pulse rectifies electrons to the conduction band in a linear section of the dI\/dV\u00a0curve, while a time-delayed weak THz pulse modulates the THz peak voltage. This approach allows waveform characterization without requiring complex retrieval algorithms31<\/a>,45<\/a> (Fig.\u2009S2e, f<\/a>). Waveform sampling at extended delays reveals THz reflections at approximately 40, 60, 80, 100, and 120\u2009ps due to sapphire windows of the vacuum chamber, cryostat, and the beam combiner, which are unavoidable in our setup (Fig.\u2009S2g<\/a>). Thickness variation of sapphire windows and the beam combiner effectively reduces the amplitude of major THz reflection to below 30% of the main pulse, causing reflections at 40 and 80\u2009ps to show distorted waveforms due to interference effects.<\/p>\n Theoretical model<\/p>\n The dynamics observed in the tunnel junction can be well represented by a time-dependent Markov process, which is described by the so called Master equation \\(\\dot{N}=M\\cdot N\\), where N is a vector containing the occupation probabilities of all considered states of the system and M is the transition matrix. Here, we considered three distinct electronic states: The electronic ground state at energy E0\u2009=\u20090\u00a0eV, where the defect is charge neutral (Q0\u2009=\u20090\u00a0e), as well as two negatively charged states (Q1,2\u2009=\u2009\u2212e), with an additional electron in either the LUMO (E1) or the LUMO+1 (E2). In addition, we model each electronic state with a single vibronic mode with 6 (14) vibrational excitation levels plus the vibrational ground state, yielding a total of 21 (45) different vibronic states for 2\u2009ML (1\u2009ML) WSe2. The different number of required vibrational levels is due to different Huang-Rhys factors between 1\u2009ML and 2\u2009ML WSe2 (see below).<\/p>\n Transition between these states can occur due to inelastic charge transfer between defect and tip or sample, or via vibrational relaxation within one electronic state. Accordingly, the transition matrix can be written as the sum of the transition rates due to sample, tip and phonon relaxation:<\/p>\n $$M={M}^{{{{\\rm{s}}}}}+{M}^{{{{\\rm{t}}}}}+{M}^{{{{\\rm{ph}}}}}$$<\/p>\n \n (1)\n <\/p>\n The transitions related to a change of charge state are given by<\/p>\n $${M}_{f{\\lambda }^{{\\prime} }i\\lambda }^{{{{\\rm{s}}}}}={T}_{f{\\lambda }^{{\\prime} }i\\lambda }({z}_{s},0)\\cdot {V}_{{\\lambda }^{{\\prime} }\\lambda }\\cdot {\\Upsilon }_{fi}$$<\/p>\n \n (2)\n <\/p>\n $${M}_{f{\\lambda }^{{\\prime} }i\\lambda }^{{{{\\rm{t}}}}}={T}_{f{\\lambda }^{{\\prime} }i\\lambda }({z}_{t},V)\\cdot {V}_{{\\lambda }^{{\\prime} }\\lambda }\\cdot {\\Upsilon }_{fi}\\cdot {\\Lambda }_{fi},$$<\/p>\n \n (3)\n <\/p>\n with indices i and f representing the initial and final electronic state, \u03bb and \\({\\lambda }^{{\\prime} }\\) the initial and final vibrational state, \\({V}_{{\\lambda }^{{\\prime} }\\lambda }\\) the Franck\u2013Condon factors, and \\({T}_{f{\\lambda }^{{\\prime} }i\\lambda }\\) the tunneling matrix element, respectively. \u03d2fi captures the effect of various multiplicities of the involved states, and \u039bfi represents their relative coupling strength to the tip.<\/p>\n The tunneling matrix element \\({T}_{f{\\lambda }^{{\\prime} }i\\lambda }(z,V)\\) describes the tunneling probability of an electron between the tip or sample and the defect state, depending on bias voltage V and distance z between the defect and lead:<\/p>\n $${T}_{f{\\lambda }^{{\\prime} }i\\lambda }(z,V)=\\left\\{\\begin{array}{ll} \\int_{-\\infty }^{V}{{{\\rm{d}}}}E\\,g(E-\\Delta {E}_{f{\\lambda }^{{\\prime} }i\\lambda })\\,{e}^{-2\\kappa z}\\quad &,{{{\\rm{if}}}}{Q}_{f}-{Q}_{i}=-e\\\\ \\int_{V}^{\\infty }{{{\\rm{d}}}}E\\,g(E+\\Delta {E}_{f{\\lambda }^{{\\prime} }i\\lambda })\\,{e}^{-2\\kappa z}\\quad \\;\\;\\; &,{{{\\rm{if}}}}{Q}_{f}-{Q}_{i}=+ e\\\\ 0\\, \\qquad \\qquad\\qquad\\qquad\\quad\\qquad\\quad &,\\; {{{\\rm{otherwise}}}} \\;\\;\\quad \\end{array}\\right.$$<\/p>\n \n (4)\n <\/p>\n with \\(\\Delta {E}_{f{\\lambda }^{{\\prime} }i\\lambda }\\) being the energy difference between the initial and final vibronic state and decay rate \\(\\kappa=\\sqrt{2{m}_{e}(\\phi -E+eV\/2)}\/\\hslash\\). The transition threshold is broadened by a Gaussian line shape g(E) and normalized to integrate to the quantum of conductance G0\u2009=\u20092e\/\u210f for z\u2009=\u20090.<\/p>\n For the Franck\u2013Condon factors \\({V}_{{\\lambda }^{{\\prime} }\\lambda }\\), we assume only a rigid shift in normal mode coordinates and negligible Duschinsky rotation46<\/a>:<\/p>\n $${V}_{{\\lambda }^{{\\prime} }\\lambda }=\\frac{\\min ({\\lambda }^{{\\prime} },\\lambda )!}{\\max ({\\lambda }^{{\\prime} },\\lambda )!}\\,{S}^{| {\\lambda }^{{\\prime} }-\\lambda | }\\,{e}^{-S}\\,{\\left[{L}_{\\min ({\\lambda }^{{\\prime} },\\lambda )}^{| {\\lambda }^{{\\prime} }-\\lambda | }(S)\\right]}^{2},$$<\/p>\n \n (5)\n <\/p>\n with Huang-Rhys factor S and Laguerre polynomials \\({L}_{n}^{\\alpha }\\).<\/p>\n The effect of the various multiplicities on the transition rates is captured by the factor \u03d2fi. Since we are only considering transitions for which either the final or initial state is a singlet, we can write \u03d2fi\u2009=\u2009mf, with mf being the multiplicity of the final state.<\/p>\n The last term \u039bfi for the tip-mediated transitions captures the relative coupling of LUMO and LUMO+1 to the tip, which can vary with lateral tip position. The additional factor \u039bfi is normalized to 1 for the tip\u2013LUMO coupling and has typical values between 0.5 and 1.0 for the tip\u2013LUMO+1 coupling, chosen to match the experimental dI\/dV spectra.<\/p>\n For the charge-neutral phonon relaxations, the transition matrix is written as a decay to the vibrational ground state with lifetime \u03c4ph:<\/p>\n $${M}_{f{\\lambda }^{{\\prime} }i\\lambda }^{{{{\\rm{ph}}}}}=\\left\\{\\begin{array}{ll}1\/{\\tau }_{{{{\\rm{ph}}}}}\\quad &{{{\\rm{,\\; if}}}}f=i,{\\lambda }^{{\\prime} }=0{{{\\rm{and}}}}\\lambda > 0\\\\ 0 \\qquad\\quad &,{{{\\rm{otherwise.}}}} \\qquad\\qquad\\;\\;\\end{array}\\right.$$<\/p>\n \n (6)\n <\/p>\n Lastly, the diagonal entries of the total transition matrix are set to Mll\u2009=\u2009\u2212\u2211k\u2260lMkl.<\/p>\n For static bias voltages, the master equation is solved for \\(\\dot{N}=0\\) to obtain the bias and z dependent equilibrium occupation Neq. The dc current is then given by the net charge transfer between defect and sample \\({I}_{{{{\\rm{dc}}}}}={\\sum }_{k}{\\left({W}^{s}N\\right)}_{k}\\), with \\({W}_{f{\\lambda }^{{\\prime} }i\\lambda }^{s}={M}_{f{\\lambda }^{{\\prime} }i\\lambda }^{s}\\cdot ({Q}_{f}-{Q}_{i})\\). In order to evaluate the evolution of the system under a transient bias voltage, the master equation is solved in discrete time steps \\({N}_{i+1}-{N}_{i}=M\\left(V({t}_{i})\\right)\\Delta t\\cdot {N}_{i}\\), starting from the equilibrium occupation at Vdc.<\/p>\n Implementation of simulation<\/p>\n The simulations are\u00a0based on the master equation\u00a0and require experimental parameters that can be obtained from static measurements of the VacSe (Fig.\u2009S1<\/a>). We extracted the energies of localized defect states E1 and E2 via the zero-phonon lines and the tip coupling factor Lf,i from the experimental dI\/dV spectrum at the set point z0. The average charge-state lifetime \u03c40, and the tunneling barrier \u03ba were obtained via a rate equation fit to the I(z) approach curve16<\/a>.<\/p>\n The vibronic broadening of the defect state resonances was approximated by a single vibrational mode with \u210f\u03a9\u2009=\u20098\u2009meV and a Huang-Rhys factor of S\u2009=\u20092.2 (2\u2009ML) for both charge states. We find that a Gaussian lineshape with 3\u2009meV bandwidth provides a good match with the experimental dI\/dV (Fig.\u2009S6b<\/a>), while being consistent with literature on VacS in WS2, observing a dominant phonon mode of 12\u2009meV and a Huang-Rhys factor of 1.218<\/a>. Considering an average 40% softening of the phonon modes in WSe2 as compared to WS247<\/a>, we choose a phonon mode of \u210f\u03a9\u2009=\u20098\u2009meV, and fit the Huang-Rhys factor to the STS spectrum. For the simulations on 1\u2009ML WSe2 we cannot directly fit the vibronic satellite peaks of the LUMO resonance. Instead, we rescale the Huang-Rhys factor from 2\u2009ML WSe2 by a factor of about two, as previously found for CS in 1\u2009ML and 2\u2009ML WS219<\/a>. Therefore, we choose S\u2009=\u20095 for VacSe in 1\u2009ML WSe2, while we assume the phonon mode energy to remain constant, \u210f\u03a9\u2009=\u20098\u2009meV19<\/a>.<\/p>\n Although lifetime broadening of a state with tens of picoseconds of lifetime would lead to a Lorentzian lineshape with a few tens of \u03bceV broadening, we justify the choice of Gaussian lineshapes by the presence of additional low-energy vibrational modes with a high Huang-Rhys factor (S\u2009\u226b\u20091), e.g., well-known interlayer phonons for multi-layer TMDs48<\/a>. These low-energy phonons impose an additional Franck\u2013Condon blockade that explains the difference between simulation and experiment at bias voltages overlapping with the LUMO resonance (\u0394V\u2009>\u20090\u00a0V). In particular, the breakdown of simulations at \u0394V\u2009\u2248\u200910\u2009mV in Fig.\u20093<\/a>c emphasizes the limitation of the 1-phonon mode treatment. A thorough 2-mode treatment would require significantly more computational resources without major benefits to our model.<\/p>\n The effective strength of the Franck\u2013Condon blockade also depends on the phonon lifetime \u03c4ph. We distinguish three regimes in the evolution of charge-state transitions between tip and defect: (I) A voltage gate for charging, (II) vibronic relaxation, and (III) charge relaxation, sketched in Fig.\u2009S3<\/a>. The THz peak voltage \\({V}_{{{{\\rm{THz}}}}}^{{{{\\rm{pk}}}}}\\) of few-100\u2009mV enables transitions from the neutral ground state to high vibronic modes of the charged state. Vibronic relaxation in a 2D semiconductor with strong electron-phonon coupling typically occurs faster than the THz gate timescale. Additionally, due to local resolution of LW-STM, combined with the fast lateral diffusion of phonon modes, the atomic-scale QD remains unaffected by the lifetime of acoustic phonons. As a result, we expect relaxation of the charged VacSe to the vibrational ground state within few-100\u2009fs, as seen for thermalization and cooling of photoexcited charge carriers37<\/a>. We conservatively assume \u03c4ph\u2009=\u20091\u2009ps for the simulations. After the source pulse and at small \u0394V\u2009\u2272\u20090 V, charge-state transitions from low vibrational modes of the charged state have only a few vibrational modes in the neutral state available, creating an asymmetry between forward and backward tunneling that reduces the impact of back tunneling. In addition, our simulation considers the LUMO and LUMO+1 multiplicity as a quartet (J\u2009=\u20093\/2) and a sextet (J\u2009=\u20095\/2) state, respectively18<\/a>, further enhancing this asymmetry (Fig.\u2009S5<\/a>).<\/p>\n The voltage drop in 2\u2009ML WSe2 is approximated by a plate capacitor model, \\({V}_{{{{\\rm{eff}}}}}={V}_{{{{\\rm{dc}}}}}\\frac{z}{z+d\/{\\epsilon }_{r}}\\), with a gap composed of vacuum and WSe2 using a dielectric constant \u03f5r\u2009=\u20096.5 and d\u2009=\u20096.5\u2009\u00c5 for WSe2 layer spacing49<\/a> (Fig.\u2009S6e<\/a>). Due to the complex tip shape and iso-potential surfaces that are not parallel to the substrate plane, the voltage drop\u00a0does not only depends on the tip height z but also on the lateral position of the tip. An accurate description of the effective bias potential can be obtained by integrating the wavefunction with the spatial dependence of the electric potential within the junction50<\/a>. In most simulations shown in this work, this effect is not crucial to reproduce the experimental data. However, for the QLW(z) measurements, this simple plate capacitor model is not sufficient to model the experiment quantitatively, see Extended Fig.\u2009S6g<\/a> (red line). Because the Franck\u2013Condon blockade depends strongly on the effective applied Vdc, it is very sensitive to changes in the voltage drop and therefore needs a more sophisticated model for the voltage drop. For this reason, we introduced a slightly modified junction model for the QLW(z) measurements. The tip position for tunneling is assumed to be at height z0, whereas the effective plane of the plate capacitor is at a larger distance from the sample zV D\u2009>\u2009z0 (Fig.\u2009S6e\u2013g<\/a>). For the 1\u2009ML WSe2 the simple plate capacitor model is used with an d\/\u03f5r\u2009=\u20090.1\u2009\u00c5.<\/p>\n