Ultrafast Coulomb blockade in an atomic-scale quantum dot

Domain: physics
Published: 2025-12-01
Editorial outcome: accept

Abstract
Controlling electron dynamics at optical clock rates is a fundamental challenge in lightwave-driven nanoelectronics and quantum technology. Here, we demonstrate ultrafast charge-state manipulation of individual selenium vacancies in monolayer and bilayer tungsten diselenide using picosecond terahertz source pulses, focused onto the junction of a scanning tunneling microscope. Using pump–probe time-domain sampling of the defect charge population, we capture atomic-scale snapshots of the transient Coulomb blockade, a hallmark of charge transport via quantized defect states. We leverage the Franck–Condon blockade, which restricts accessible vibronic transitions and promotes unidirectional charge transport, to effectively mitigate back tunneling to the tip electrode. Our master equation approach models the non-reciprocal tunneling process due to vibrations and angular momentum multiplicities, accurately reproducing the time-dependent tunneling current across different coupling regimes. Capturing and controlling ultrafast charge dynamics in low-dimensional materials at the atomic scale opens frontiers in lightwave-driven nanoscalescience and technology.

Introduction
Artificial atom qubits, such as quantum dots (QDs) and color centers, defined atomic interfaces are two-dimensional (2D) materials that have attracted significant interest in light of their applications in are chemically versatile, tunable via external fields, and they can host quantum sensing, communication, and photonic quantum atomic-scale QDs9. Several spin-bearing, optically active defect emit- computing1–3. In particular, the inherently strong confinement for ters have been identified in transition metal dichalcogenides solid-state atomic-scale defects imposes a large characteristic energy (TMDs)10,11 and hexagonal boron nitride12,13, including substitutional scale that leads to exceptionally long coherence times, and operability defects and vacancies. However, achieving device-scale electron beyond cryogenic environments4. Spin-selective optical decay path- transport through single defects in 2D materials has so far remained ways enable high-fidelity spin initialization and optical readout5. While elusive.

The tunneling contact in scanning probe microscopy offers a optical access is advantageous for addressing individual QDs at THz to means of studying quantized electron transport through single defects PHz frequencies6, electrical interfaces offer crucial benefits for scalable and exploring the Coulomb blockade regime at the atomic scale. on-chip integration. The primary engineering challenge lies in devel- Steady-state Coulomb blockade of adatoms, molecules, and defects oping suitable leads for atomic interfaces. Although extreme ultra- has previously been observed in scanning tunneling microscopy violet lithography is making progress towards atomic contacts, (STM)14–16. Dangling bonds of an H-terminated Si surface have been scanning probe hydrogen resist lithography7,8 has already demon- demonstrated to exhibit quantized transport8,17 and Coulomb strated proof-of-concept devices. Promising candidates for well- blockade8. Similarly, chalcogen vacancy defects in TMDs introduce coupling between layers and across 2D heterostructures results in discrete defect states deep in the band gap of the host material, with long-lived excitations. As a result, the regime of QLW > 1e/pulse is not narrow defect resonances followed by vibronic sidebands18,19. These always accessible, because the large Coulomb energy associated with defect states can be probed individually via scanning tunneling spec- double charging of the QD inhibits multiple electron transfers within troscopy (STS), revealing the influence of spin-orbit coupling18 and the short pulse duration. However, by using two time-delayed THz layer-dependent binding energy16,20. To operate and control single- pulses, it becomes possible to probe e−–e− correlations, such as the electron dynamics in these atomic QDs at their natural time scale, Coulomb blockade in the time domain.

Previous studies could show that the average charge-state lifetime τ0 can be estimated via the maximum tunneling rate at saturated current through a localized state15,32. Specifically, when the tip–sample coupling ΓT greatly exceeds the substrate coupling ΓS, the relationship Γ−1 = Γ−1 + Γ−1 ≈ τ0 holds, which can be achieved by approaching the tot S T tip. Additional WSe2 layers effectively decouple VacSe and substrate, allowing control over τ0 across several orders of magnitude. In a previous study, we reported average charge-state lifetimes of single VacSe in the top layer of 1–4 ML WSe2, ranging from 1 ps to several ns16. Owing to the spatial distribution of the defect orbital, top and bottom VacSe in the same layer exhibits similar coupling strength to the substrate and thus comparable charge-state lifetime16. This study lays the basis for investigating ultrafast charge dynamics in the time domain using LW-STM. As a starting point, we consider a neutral VacSe in 1 ML and 2 ML WSe2 characterized by the LUMO and LUMO+1 states (Fig. 2a), residing in the unoccupied spectrum of the WSe2 bandgap. STM topography and orbital images reveal the three-fold symmetric orbital shape of the defect state (Fig. 2b)18.

Quenching of lightwave-driven orbital imaging
To maximize the time resolution of THz pulses, the majority of LW-STM/STS experiments have so far been performed in a configuration where Vdc ≪ VTHz and at or close to Vdc = 0 V to minimize dc background currents. Figure 2c shows the rectified charge per THz waveform (transient) in units e/tr as a function of peak THz voltage, which, in first order, reproduces the I(V) characteristic of the tunnel junction28. The onset of the two rising edges in the QLW(VpkTHz) curve (black arrows) aligns with the expected spectral locations of the lowest unoccupied molecular/defect orbital (LUMO) and LUMO+1. Additionally, the LW-STM map of the LUMO [Fig. 2d (i)] closely resembles the dc STM measurement. However, the situation is strikingly different for LW-STS of VacSe on 2 ML WSe2 (Fig. 2c turquoise). In this case, we are unable to resolve the defect states through LW-STM [Fig. 2d (iii) and (iv)], even though their STS spectrum matches closely the spectrum of VacSe on 1 ML WSe2. Instead of observing an enhancement in the LW-driven rectification at the defect position, which would typically arise from additional tunneling channels, we see a reduction of LW-driven tunneling at the defect center. This observation suggests that the VacSe on 2 ML WSe2 enters a different lifetime regime, affecting the ultrafast tunneling process. The sub-nm point-like reduction of LW rectification to the conduction band in configurations (ii) and (iv) differs from the defect orbital and appears as a consequence of local band bending due to the charged VacSe33,34.

We attribute the loss of in-gap defect state contrast in 2 ML to an enhanced charge-state lifetime, related to the electronic decoupling of VacSe via the additional WSe2 layer16. The average charge-state lifetime of different VacSe/2 ML obtained from dc current saturation measurements ranges from τ02ML = 50 to 86 ps (Fig. S1), significantly longer than the sub-1 ps source voltage duration. Consequently, multiple tunneling events per pulse are highly unlikely. Subsequently to the injection by the source pulse, the e− can relax through two possible pathways: (i) forward tunneling to the substrate, contributing to a measurable current, or (ii) backward tunneling to the tip. Due to the small duty cycle in LW-STM, intermediate to strong charge injection from the tip is required to facilitate tunneling within the short duration of the THz peak field. Hence, back tunneling is the more probable pathway if the charge-state lifetime exceeds a few picoseconds, quenching charge rectification. As a result, a transiently charged VacSe in 2 ML is more likely to relax by discharging via the STM tip (back tunneling) instead of the substrate (forward tunneling), resulting in zero net charge rectification. In contrast, the smaller lifetime of VacSe in 1 ML (τ01ML < 3 ps) enables notable charge rectification from tip to substrate during the transient THz pulse, which explains the orbital contrast observed in the LW-STM images for VacSe in 1 ML WSe2. For THz peak voltages reaching the conduction band, femtosecond carrier relaxation and strong delocalization of continuum states enable sub-cycle charge rectification during the THz gate, irrespective of the layer thickness, yielding an unperturbed onset of conduction band states (gray curves in Fig. 2c).

Back tunneling suppression via Franck–Condon blockade
Next, we study the effect of the dc bias on THz charge rectification for a VacSe in 2 ML (Fig. 3a). To ensure comparability across different VacSe, we define a relative voltage offset from the LUMO ΔVz = Vdc − VLUMO(z) to account for local shifts of the relative LUMO position with tip height, lateral tip position, and sample heterogeneity (Methods). The first two factors arise from the well-known voltage drop across a double-barrier tunneling junction35. In contrast to Fig. 2d (iii), where we set Vdc = 0 V, we start to observe orbital contrast in LW-STM as Vdc nears the LUMO resonance, as seen in Fig. 3b. This effect becomes more pronounced as ΔVz approaches zero. We attribute the reestablished LW-STM contrast to a suppression of back tunneling. This suppression can be explained by the so-called Franck–Condon blockade36, which is more significant than the (slowly-varying) bias-dependent change of the tunneling barrier. Note that the orbital shape of VacSe in Fig. 3b appears saturated, unlike in the monolayer case in Fig. 2d (i). Due to the required ΔVz ≈ 0 V, dc electrons block the THz electrons at some tip positions. The same orbital shape as for the 1 ML case can be recovered for specific Vdc and tip heights, as shown in Fig. 4b. When Vdc exceeds the onset of the VacSe resonance (ΔVz > 0 V), significant dc currents occur, eventually leading to the saturation of Idc and causing charge fluctuations that substantially increase the noise in the THz current (Fig. 3c).

In the Franck–Condon blockade regime, changes in the phonon population of a defect can introduce an asymmetry between forward and backward tunneling. In Fig. 3d, e, the principle of the Franck–Condon blockade is schematically illustrated. Within the Franck–Condon approximation, each charge-state of the defect is represented by a harmonic potential with a single vibrational mode of characteristic energy ℏΩ. Due to fast vibrational cooling in 2D semiconductors37 (curved black arrow), charge-state transitions primarily occur from the vibrational ground state (dark red), centered around the relaxed nuclei positions. The displacement of the potential energy surface between different charge states increases with the electron-phonon coupling strength. As a result, electron transfer (vertical blue arrows) is generally accompanied by vibrational excitation. However, when the dc bias is set close (within the order of the vibrational energy) to the LUMO orbital energy (Fig. 3d), the overlap with energetically accessible low-lying vibrational states of the neutral defect is greatly reduced38. In contrast, when the dc bias is much lower (Fig. 3e), all vibrationally excited states remain accessible. Consequently, back tunneling can be significantly suppressed by raising the electrochemical potential of the tip.

We model the transient charging and discharging of the atomic defect by a time-dependent Markov process, described by the so-called Master equation, as detailed in the Methods section. This model treats the defect by three electronic states (one ground state and two charged states, LUMO and LUMO+1), as well as a single vibrational mode (ℏΩ ≈ 8 meV) with multiple quanta, and a Huang-Rhys factor of S ≈ 2.2 for a VacSe in 2 ML WSe2. These parameters are extracted from experimental STS spectra and the average charge-state lifetime τ0 from STM approach curves (Fig. S6a, b). The transition matrix is expressed as the sum of the transition rates due to sample, tip, and phonon relaxation. To calculate the rectified charge, we model the actual waveform of the THz transient (Fig. S6c). The simulations, shown in Fig. 3c (red), qualitatively and quantitatively reproduce the steep increase of the THz current as a function of dc bias within tens of meV of the defect resonance. Importantly, the strong suppression of back tunneling at elevated bias is attributed to a reduction in the energetically allowed vibronic transitions into the neutral state. The Franck–Condon blockade acts as an effective backflow valve, favoring unidirectional charge transport when charge-state lifetimes exceed the pulse duration, which is crucial for achieving charge rectification in LW-STM.

Optimal tip distance between weak and strong coupling
The coupling of the atomic QD to the tip contact not only sensitively depends on the tip bias, but also on the tip–sample distance (Fig. 4a). Figure 4b shows a series of LW-STM images illustrating the evolution of the orbital contrast at different tip heights at Vdc = 0.7 V, corresponding to ΔVz = 0 at z0. In contrast to Fig. 2c the orbital images are saturated as a consequence of a transient population that produces a convolution of defect orbital and a local variation in Franck–Condon blockade. The black curve in Fig. 4c presents a QLW approach curve obtained on a VacSe orbital lobe. Both spatial (Fig. 4b) and height (Fig. 4c) dependence reveal an initial increase, followed by a decrease of the LW signal as the tip approaches the defect. Initially, reducing the tip height compensates for the short duty cycle in LW-STM and enhances the signal. However, at very close distances, stronger back tunneling suppresses the rectified current, resulting in a maximum rectification near z − z0 = −2.5 Å. Two factors contribute to this effect: Due to the voltage drop across the WSe2 layers, smaller tip heights (or closer lateral tip position) shift the LUMO resonance to higher voltages, thereby reducing ΔVz and thus enhancing back tunneling, as previously discussed. Additionally, the increasing tunneling rate ΓT(z) becomes much larger than the defect–graphene tunneling rate, ΓS. This promotes back tunneling, dominating over the Franck–Condon blockade. Our simulations (Fig. 4c, red) reproduce the shape and magnitude of the experimental QLW(z) curve, accounting for both the z-dependent voltage drop and the increase in ΓT(z). While quenching of excited states via electron transfer to the tip was previously reported for photo-excited molecules39,40, here we directly map this process with ultrafast time resolution.

Time-domain detection of the ultrafast Coulomb blockade
Finding optimal parameters for ΔVz and z enable time-domain mapping of the transient Coulomb blockade close to its average charge-state lifetime. Here, we use a sequence of two time-delayed THz transients in combination with Vdc, such that Vdc + VpkTHz accesses the localized in-gap defect state, but remains below the onset of the conduction band (Fig. 5a, d). In Fig. 5b, e, we show the transient Coulomb blockade at VacSe in 1 and 2 ML WSe2. In these measurements, QLW(Δt) describes the rectified charge by the probe transient only (kHz amplitude modulation), while the pump transients interact at full repetition rate (41 MHz). Near the overlap Δt = 0 ps, the signal is dominated by interference between the THz waveforms, such that we restrict our analysis to time delays Δt > 1.5 ps. For time delays much longer than the average charge-state lifetime τ0, the transient signal recovers to the equilibrium value Qeq, corresponding to unperturbed THz rectification by the probe transient. For orbital images acquired at different time delays, we plot the relative change of rectified charge ΔQLW = QLW − Qeq, obtained by subtraction of a reference map at a large delay (Δt ≫ τ0). The differential images clearly reveal the orbital structure of VacSe with negative amplitude, providing a direct and dynamic visualization of a Coulomb blockade at the atomic scale.

We model the lightwave-driven pump-probe (LW-PP) point spectrum by a single exponential fit (dashed black line) that allows to extract the relaxation dynamics of the VacSe QD. This effective charge-state lifetime τeff combines the interplay of intrinsic charge relaxation to the substrate and ultrafast back tunneling to the tip. To disentangle these contributions, we record LW-PP spectra at increasing Vdc. We determine τeff via exponential fitting of the spectra at time delays between 1.5 ps and 12 ps, ensuring robust convergence of fit parameters. The error bar includes disturbances related to THz trailing oscillations at high Vdc. τeff strongly increases as ΔVz → 0 and peaks near ΔVz = 0 V, where it reaches the intrinsic charge-state lifetime τ0 (Fig. 5c). At ΔVz ≈ 0 V, back tunneling is strongly suppressed and the intrinsic relaxation channel to the substrate becomes the dominant mechanism for charge relaxation. Notably, the effective charge-state lifetime is rapidly quenched for negative ΔVz. This highlights the critical role of the Franck–Condon blockade for this experiment. At ΔVz < −50 mV the time resolution and transient signal amplitude become insufficient to resolve the Coulomb blockade. Our master equation model enables simulations of charge rectification at various time delays using a model of our experimental waveform (Methods). Exponential fitting of simulated LW-PP spectra with similar boundary conditions as in the experiment clearly reveals the trend of increasing effective charge-state lifetime near ΔVz = 0 V (Fig. 5c). The onset of quasi-static charging at dc biases above the LUMO resonance, as shown in Fig. 3c, precludes ultrafast measurements due to dominant stochastic tunneling.

We performed the same series of measurements at a VacSe in 2 ML WSe2 with an intrinsic charge-state lifetime τ0 = 86 ps. As the charge-state lifetime exceeds the time window of trailing reflections from windows in our beam path, which can amount up to 30% field amplitude (Fig. S2), QLW(Δt) deviates from a simple exponential decay, instead displaying dips indicative of interference with these reflections. Since ΔVz values close to 0 V are necessary to suppress back tunneling effectively, the impact of trailing reflections on the measurements can only be minimized, though not entirely eliminated. At elevated Vdc close to the onset of the defect state, the amplitude of THz pulse reflections is sufficiently high to interact with the VacSe states, causing additional excitations and experimental replica at various time delays. Nevertheless, an ultrafast Coulomb blockade lasting over few-10 ps is clearly revealed by the LW-PP spectrum and ultrafast orbital images in Fig. 5e. For simple and robust extraction of an effective charge-state lifetime, we determine the equilibrium signal at Δt > 250 ps and fit a single exponential model only to the unperturbed data from 1.5 to 15 ps. Simulated LW-PP spectra on the basis of the Franck–Condon blockade and the full transient waveform with multiple reflections quantitatively model the experiment, including the effect of THz pulse reflections (Methods, Fig. S2). Limiting the fit interval to Δt < 15 ps causes large fit errors for τeff ≳ 20 ps, however, it ensures robustness of the fit in experiment and simulation. Again, the sharp decrease of τeff with more negative ΔVz emphasizes the increasing relevance of back tunneling for VacSe in 2 ML with extended τ0. As ΔVz strongly depends on the lateral and vertical position of the tip, distinct ΔVz regimes can coexist within a single LW-STM image. In the defect periphery, where ΔVz > 0, large dc currents permanently charge the defect, effectively quenching QLW. At the circular transition region, where ΔVz ≈ 0, charge fluctuations dominate, leading to significant noise in QLW (Fig. S5c).

The lack of a WSe2 moiré structure in our sample reduces mechanical strain and facilitates fast delocalization of phononic excitations. At time scales < 5 ps, signatures of coherent motion of the WSe2 layer interfere with THz field-induced modulations and remain indistinguishable in the experiment. In consequence, we do not resolve a coherent motion of the WSe2 layer or energy modulation of the defect state driven by inelastic excitations of the VacSe states as reported by Roelcke et al.28, which would also modulate the Franck–Condon blockade significantly. In contrast to fast charge relaxation on a metallic substrate, finite charge-state lifetimes and the Coulomb blockade in our samples impose a stochastically uncorrelated perturbation to any coherent signals since the charge relaxation is not phase-locked. Furthermore, the decoherence of lattice excitations is amplified by charge transfers from pulse replicas that become relevant at elevated Vdc.

In conclusion, using THz pump–THz probe time-domain sampling, we capture real-space snapshots of the transient Coulomb blockade. By fine-tuning the transient THz waveform, dc bias offset, and tip height, we selectively populate specific defect orbitals of single selenium vacancies in monolayer and bilayer WSe2 and track the charge transfer into the substrate. A rate equation model accurately describes the time-dependent tunneling process across the different coupling regimes, applicable to other materials and quantum systems. Our experiments overcome challenges posed by back tunneling, revealing a Coulomb blockade that persists for two orders of magnitude longer than the duration of THz voltage pulses. Quenching of the LW-STM signal due to population depletion by the tip electrode is mitigated by breaking the symmetry between charge transfer from and to the tip via the Franck–Condon blockade, acting as a charge backflow valve. The Franck–Condon blockade drastically reduces the charge-transfer rate from the defect state to the tip by limiting the number of available vibrational transitions. Precise control of the transient bias and coupling strength between the STM tip and the defect state enables net charge rectification, essential for LW-STM measurements. This work complements recent advances in observing defect resonances29, coherent lattice dynamics28, charge density waves41, and bidirectional subcycle tunneling currents42 with LW-STM, establishing a pathway for probing ultrafast charge dynamics in low-dimensional materials and localized quantum states. We envision that tunneling asymmetries induced by vibrational transitions, spin and orbital momentum multiplicities will enable lightwave-driven charge control in complex materials, pushing the frontiers of ultrafast nanoscale electronics.

Results
Coulomb blockade at a single Selenium vacancy
The STM geometry enables the resonant exchange of electrons (e−) in the single-atom VacSe QD between tip (source electrode) and sample (drain) through sequential tunneling processes (Fig. 1a, b). A finite tunneling current occurs if the e− exchange proceeds in series, i.e., from the source (tip) to the QD to the drain (substrate), or vice versa. Thereby, the coupling strength to the source and drain can be precisely controlled via the tip–sample distance and TMD layer thickness, respectively. In standard STM/STS, static lead potentials enable measuring energy-selective, stochastic tunneling processes and average QD populations. Ultrafast pump-probe schemes in LW-STM allow us to access also the charge dynamics of the atomic QD. In the near-field of the tunneling junction, strong-field THz source pulses act as ultrafast bias transients that add to any static dc bias, as illustrated schematically in the equivalent circuit diagram in Fig. 1c. We preserve state-selectivity of the tunneling process by tailoring the amplitude and carrier-envelope phase of the THz transient to a single-cycle, unipolar pulse resonant with an electronic defect state27,30,31.

The key observable in LW-STM is the rectified charge QLW transferred between tip and substrate per lightwave (THz) transient. Similar to the current in conventional STM, QLW can be efficiently increased by reducing the tip-sample distance z. In a QD geometry, a rectified charge requires two consecutive tunneling events from one lead to another via the discrete QD state. Here, two regimes can be distinguished: QLW ≤ 1e/pulse and QLW > 1e/pulse. In the latter regime, multiple e− tunnel per source pulse, which is only possible if the state that facilitates tunneling is short-lived, indicating strong QD–substrate coupling, ΓS, where the charge-state lifetime (τ0) is much shorter than the source pulse duration (τLW), τ0 = ΓS−1 ≪ τLW. The coupling of the tip electrode can be tuned continuously from the weak-injection regime, (ΓT−1 ≫ τLW) to the strong-injection regime. In van der Waals materials, and in particular for localized in-gap states (orbitals), the weak coupling between layers and across 2D heterostructures results in long-lived excitations. As a result, the regime of QLW > 1e/pulse is not always accessible, because the large Coulomb energy associated with double charging of the QD inhibits multiple electron transfers within the short pulse duration. However, by using two time-delayed THz pulses, it becomes possible to probe e−–e− correlations, such as the Coulomb blockade in the time domain.

Methods
Sample preparation
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,43. Key advantages of epitaxially grown samples on quasi-free-standing graphene (QFEG) are the electronic homogeneity (no moiré), a Fermi energy pinned in the center of the large WSe2 band gap, and weak hybridization with the substrate16. 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.

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 2s at 0.12 kV acceleration bias, 0.16 kV discharge voltage, and a sub-μA sputter current16. Multiple VacSe in 1–2 ML WSe2 were investigated for this study, emphasizing the robustness and repeatability of the technique. Fig. S1 lists STM topography, dI/dV spectra, and current saturation curves of the relevant VacSe.

STM height setpoint and average charge-state lifetime
For a reproducible definition of the tip height z0, the STM feedback loop was opened on pristine WSe2 at a setpoint current of 100 pA with applied bias voltage of 1.5 V (1 ML) or 1.2 V (2 ML), corresponding to an energy 200 meV above the conduction band onset. We estimate an absolute tip–sample distance at z0 of (7±2) Å. The average charge-state lifetime (τ0) of VacSe states scales with the decoupling layers of the sample and is determined via the current saturation in STM approach curves15,16,32. For simplicity, we consider only the first unoccupied defect state (LUMO) of VacSe for the determination of τ0. Reduced localization due to a lower binding energy slightly reduces the average charge-state lifetime of the LUMO+1 defect state16. For different VacSe, τ0 slightly vary due to differences in their dielectric environment, however, the energy spectra and orbital shape remain comparable within a few-10 meV shifts (Fig. S1).

Voltage drop and definition of ΔVz
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 setpoint z0. The relative voltage ΔVz = Vdc − VLUMO(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 ΔVz < 0 the dc current is negligible, ΔVz > 0 V generates a notable dc current and thus increases the VacSe population statically in the experiment. In orbital imaging, we define ΔVz at a single point on the orbital lobe of VacSe, avoiding complications of the lateral dependence. We determine VLUMO(z) directly from a dI/dV measurement at z0 or via extrapolation of VLUMO(z0) on the basis of measurements performed in ref. 16, Fig. 2. Owing to the high sensitivity towards thermal drifts on the sub-Å scale, we estimate an accuracy δΔVz ≈ 10 mV.

LW-STM setup and near-field waveform detection
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 μJ. The system achieves peak THz voltages at the STM tip up to 0.5, 1, and 2 V at 41, 20, and 10 MHz 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. Pump and probe pulses are generated independently and are recombined 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 K base temperature and 10−10 mbar pressure.

To discriminate lightwave-driven currents, we mechanically modulate the THz path at an intermediate focus using frequencies between 800 and 1000 Hz that are compatible with the 109 T impedance 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.

The THz peak voltage is calibrated against Vdc using a recently developed approach based on the strongly nonlinear onset of the WSe2 conduction band31 measured in a clean pristine region of the sample (Fig. S2a–c). On the basis of this amplitude calibration, we determine the relative time origin (Δt = 0) and temporal resolution of the experiment (380 fs) by cross correlation of equal THz fields where only the sum field rectifies charges (Fig. S2d). In a second step, we measure the transient near-field waveform at the STM tip using THz cross correlation (THz-CC)29,31. 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 curve, while a time-delayed weak THz pulse modulates the THz peak voltage. This approach allows waveform characterization without requiring complex retrieval algorithms31,45 (Fig. S2e, f). Waveform sampling at extended delays reveals THz reflections at approximately 40, 60, 80, 100, and 120 ps due to sapphire windows of the vacuum chamber, cryostat, and the beam combiner, which are unavoidable in our setup (Fig. S2g). 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 ps to show distorted waveforms due to interference effects.

Theoretical model
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 Ṅ = M · 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 = 0 eV, where the defect is charge neutral (Q = 0e), as well as two negatively charged states (Q = −e), with an additional electron in either the LUMO (E1) or 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 ML (1 ML) WSe2. The different number of required vibrational levels is due to different Huang-Rhys factors between 1 ML and 2 ML WSe2 (see below).

Transitions 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:

M = Ms + Mt + Mph (1)

The transitions related to a change of charge state are given by
Msfλ0iλ = Tfλ0iλ (zs, 0) · Vλ0λ · ϒfi (2)

Mtfλ0iλ = Tfλ0iλ (zt, V) · Vλ0λ · ϒfi · Λfi, (3)

with indices i and f representing the initial and final electronic state, λ and λ0 the initial and final vibrational state, Vλ0λ the Franck–Condon factors, and Tfλ0iλ the tunneling matrix element, respectively. ϒfi captures the effect of various multiplicities of the involved states, and Λfi represents their relative coupling strength to the tip.

The tunneling matrix element Tfλ0iλ (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:

Tfλ0iλ (z, V) = 
{ R · V 1 dE g(E − ΔEfλ0iλ ) e−2κz, if Qf − Qi = −e
{ V 1 dE g(E + ΔEfλ0iλ ) e−2κz, if Qf − Qi = +e
0, otherwise (4)

with ΔEfλ0iλ being the energy difference between the initial and final vibronic state and decay rate κ = 2m(ϕ − E + eV/2)/ℏ. The transition threshold is broadened by a Gaussian line shape g(E) and normalized to integrate to the quantum of conductance G0 = 2e/ℏ for z = 0. For the Franck–Condon factors Vλ0λ, we assume only a rigid shift in normal mode coordinates and negligible Duschinsky rotation46:

Vλ0λ = m min(λ, λ0) ! S |λ0−λ|/2 e−S L |λ0−λ| (S) (5)

with Huang-Rhys factor S and Laguerre polynomials Lαn.

The effect of the various multiplicities on the transition rates is captured by the factor ϒfi. Since we are only considering transitions for which either the final or initial state is a singlet, we can write ϒfi = mf, with mf being the multiplicity of the final state.

The last term Λfi 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 Λfi is normalized to 1 for the tip–LUMO coupling and has typical values between 0.5 and 1.0 for the tip–LUMO+1 coupling, chosen to match the experimental dI/dV spectra.

For the charge-neutral phonon relaxations, the transition matrix is written as a decay to the vibrational ground state with lifetime τph:

Mphfλ0iλ = 1/τph, if f = i, λ0 = 0 and λ > 0 (6)

Lastly, the diagonal entries of the total transition matrix are set to Mll = −∑k≠l Mkl.

For static bias voltages, the master equation is solved for Ṅ = 0 to obtain the bias and z dependent equilibrium occupation Neq. The dc current is given by the net charge transfer between defect and sample Idc = ∑ Wsfλ0iλ · (Qf − Qi), with Wsfλ0iλ = Msfλ0iλ · (Qf − Qi). In order to evaluate the evolution of the system under a transient bias voltage, the master equation is solved in discrete time steps Ni+1 − Ni = M(V(t)) · Δt · Ni, starting from the equilibrium occupation at Vdc.

Implementation of simulation
The simulations are based on the master equation and require experimental parameters that can be obtained from static measurements of the VacSe (Fig. S1). We extracted the energies of localized defect states E1 and E2 via the zero-phonon lines and the tip coupling factor L from the experimental dI/dV spectrum at the setpoint z0. The average charge-state lifetime τ0, and the tunneling barrier κ were obtained via a rate equation fit to the I(z) approach curve16.

The vibronic broadening of the defect state resonances was approximated by a single vibrational mode with ℏΩ = 8 meV and a Huang-Rhys factor of S = 2.2 (2 ML) for both charge states. We find that a Gaussian lineshape with 3 meV bandwidth provides a good match with the experimental dI/dV (Fig. S6b), while being consistent with the Huang-Rhys factor of S = 2.2 for VacSe in 2 ML WSe2. Considering the softening of the phonon modes in WSe2 as compared to WS247, we choose a phonon mode of ℏΩ = 8 meV, and fit the Huang-Rhys factor to the STS spectrum. For the simulations on 1 ML WSe2 we cannot directly fit the vibronic satellite peaks of the LUMO resonance. Instead, we rescale the Huang-Rhys factor from 2 ML WSe2 by a factor of about two, as previously found for C in 1 ML and 2 ML WS219. Therefore, we choose S = 5 for VacSe in 1 ML WSe2, while we assume the phonon mode energy to remain constant, ℏΩ = 8 meV19.

Although lifetime broadening of a state with tens of picoseconds of lifetime would lead to a Lorentzian lineshape with a few tens of μeV broadening, we justify the choice of Gaussian lineshapes by the presence of additional low-energy vibrational modes with a high Huang-Rhys factor (S ≫ 1), e.g., well-known interlayer phonons for multi-layer TMDs48. These low-energy phonons impose an additional Franck–Condon blockade that explains the difference between simulation and experiment at bias voltages overlapping with the LUMO resonance (ΔVz > 0 V). In particular, the breakdown of simulations at ΔVz ≈ 10 mV in Fig. 3c emphasizes the limitation of the 1-phonon mode treatment. A thorough 2-mode treatment would require significantly more parameters and is beyond the scope of the current model.

The Franck–Condon blockade depends on the phonon lifetime τph. 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. S3. The THz peak voltage VpkTHz of few-100 mV 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 fs, as seen for thermalization and cooling of photo-excited charge carriers37. We conservatively assume τph = 1 ps for the simulations. After the source pulse and at small ΔVz ≲ 0 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 of back tunneling. In addition, we consider the LUMO and LUMO+1 multiplicity as a quartet (J=3/2) and a sextet (J=5/2) state, respectively18, further enhancing this asymmetry (Fig. S5).

The voltage drop in 2 ML WSe2 is approximated by a plate capacitor model, V = Vdc · z/d, with a gap composed of vacuum and WSe2 using a dielectric constant ϵr = 6.5 and d = 6.5 Å for WSe2 layer spacing49 (Fig. S6e). Due to the complex tip shape and iso-potential surfaces that are not parallel to the substrate plane, the voltage drop does 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. 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. S6g (red line). Because the Franck–Condon 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 z > z0 (Fig. S6e–g). For the 1 ML WSe2 the simple plate capacitor model is used with d/ϵr = 0.1 Å. To match the simulations of the VacSe in 1 ML WSe2 in Fig. 5 quantitatively, we assume a 0.4 Å uncertainty on the experimental z0 that may arise from small lateral variations in the tip position between the I(z) spectroscopy and THz pump-probe measurements.

The simulation includes a model of the extended transient waveform to reproduce the experimental time-domain signal. Since pump and probe are indistinguishable in the experiment and equally interact with the VacSe state, the Coulomb blockade causes a transient signal that is symmetric at Δt = 0 ps. This significantly impacts measurements at Δt > 15 ps, where transient signals from the main pulse and reflections overlap, causing a non-trivial curvature. Thus, we restrict the data interval for exponential fitting between 1.5 and 15 ps. As another consequence of reflections, extended coherent signals are quenched due to destructive interference from multiple excitations.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Allerbeck et al. report the lightwave-driven manipulation of a charge state of a single atomic vacancy in mono- and bilayer\nWSe2. Building on a Franck-Condon suppression of back-tunnelling and lightwave-driven STM, they develop a new method\nthat allows tracking the occupation dynamics on atomic scales and with ps to fs temporal resolution. By probing the vacancy\nfaster than a characteristic lifetime of a charged state, they directly observe the Coulomb blockade unfolding in time. The\nauthors build an effective dynamical model to describe the role of tip-vacancy coupling strength in the Coulomb blockade\nmeasurements. The experiments are performed at the very high technical level, and the data supports the main claims of the\nmanuscript.\nBy merging the atomic-femtosecond resolution of the lightwave-driven microscopy with the direct observation of Coulomb\nblockade, this manuscript provides pivotal insights to the long-standing issue of directly measuring the charge state\ndynamics in single molecules and atomic-scale building blocks of modern nanotechnology. Mapping the dynamic effects of\ncoupling strength provides unprecedented information for designing lightwave electronics. This manuscript has the potential\nto become a journal highlight, and I am happy to recommend it for publication in Nature Communications after a minor\nrevision.\nI would ask the authors to address the following issues:\n1. One of the key novelties of the manuscript is undoubtedly the measurements of the single defect charging suppression\nwith an ultrafast temporal resolution demonstrated in Fig. 5. I find the data very convincing, yet I have minor reservations\nabout the interpretation:\n1.1. For the monolayer case, I am not yet convinced that the fit in Fig. 5b (bottom panel) starting as early as 1.5 ps is justified.\nWith DC offset bias almost aligned with the LUMO, it is the fine structure of the THz transients with all their trailing\noscillations that becomes important. Indeed, such structure is visible at least up to 2 ps in Fig. S2f and can be resolved even\nbeyond as oscillations exceeding the apparent noise level in Fig. S2g. I believe that (i) full cross-correlation measurements\nbetween pump and probe pulses at the ML vacancy orbital could be shown and (ii) the effect of the trailing\noscillations quantified.\n1.2. What is the effect of stochastic tunnelling on the ultrafast Coulomb blockade measurements when V_DC approaches\nV_LUMO (e.g. in Fig. 5c)? In Fig. S1, the authors clearly show that at z = z_0 – 3.5A and V_DC ~ V_LUMO + 0.2V the DC\ncurrent shows Coulomb blockade signatures in a purely static experiment. How does the\nstochastic occupation compare to\nthe lightwave-driven one? I believe that Fig. 5c would benefit from showing a DC current curve as a function of V_DC\ndetuning.\n1.3. Semi-log scale in exponential decay curves would be more convincing and would allow the reader better estimate the\nnoise level and fit reliability.\n2. Is it possible to image the vacancy orbital in 2ML with LW-STM (e.g. in Fig. 2d (iii)) by increasing the integration time? In\nFig. 2c, with a larger integration time, the turquoise curve is above the grey one, suggesting that the orbital should be visible\nas compared to the surroundings. If yes (e.g. a linecut across with orbital with 100x integration time), could this estimate be\nused to quantify the back-tunnelling rate?\n3. The authors discuss possible vacancy level shifts due to coherent lattice dynamics introduced by the THz pump pulse,\ncomparing their work to Ref. 28. Yet, the present data shows no out-of-plane coherent lattice motion in either the pristine ML\nor BL, in contrast to the drum mode in Ref. 28. What kind of lattice dynamics do the authors imply then? In-plane? Can the\nauthors rule out the effects of coherent acoustic phonon wave packets observed by Sheng et al. in PRL 129, 043001\n(2022)? I would also take the liberty to notice, that Ref. 28 reveals dynamic shifts of the vacancy states even in the absence\nof charging, in contrast to what is written in the present manuscript.\n4. I encourage the authors to compare their technique with other recent works in the field of lightwave-driven microscopy. In\nparticular, I would like to bring the authors’ attention to the recent paper by Siday et al., Nature 629, 329-334 (2024) that\nreports the subcycle tunnelling currents within the THz pulse itself, and is thus sensitive to the non-rectified, including the\nback-tunnelling, charges (e.g. Fig. 5 in the Siday et al.).\n5. The authors very convincingly highlight the importance of and the hopes related to the defect-based functional properties,\nso it would be great to see what the field can learn from this manuscript. Is the regime of strong coupling to the leads\naccessible in e.g. vertical heterostructures, THz or optical photonic resonators? The spatial probing of the ultrafast Coulomb\nblockade (Figs. 4, 5) is very impressive: could it provide us with some ideas for the device design or be used in device\nsimulations? Could this technique be applied for molecules? I would appreciate a broader outlook.\nManuscript presentation suggestions:\n6. Is the use of the picocavity concept really relevant? The authors mention it only once and briefly, in the abstract, without\nmuch context.\n7. In Fig. 4, the authors seem to use the term LW-STS to describe the z-dependence of the signal. While it z-spectroscopy is\nalso a kind of spectroscopy, it was confusing for me as a reader.", "reviewer_2": "The work by Allerbeck et al. demonstrates picosecond-resolved control and detection of single-electron charge dynamics in\nan atomic-scale quantum dot (VacSe defects in WSe₂) using lightwave-driven STM (LW-STM) and a THz pump–THz probe\nconfiguration. The key finding is the direct observation of ultrafast Coulomb blockade and its modulation via the Franck–\nCondon (FC) blockade. I think this is a significant achievement in the field of atomic quantum transport.\nEven though charge state lifetimes have been studied extensively using conventional STM (see cited literature), the\nscientific results enter a new regime of real-time observation of picosecond charge state lifetimes. The main novelty in this\nwork is the ability to observe the transient Coulomb blockade dynamics in real time and down to a few picoseconds. In\naddition, the authors demonstrate active control of tunneling directionality by inducing an asymmetry in forward-backward\ntunneling by the Franck–Condon Blockade, which is a new achievement, despite the fact that the underlying physics is well-\nestablished. The paper is well written, the data and figures are presented in good quality and the methodology and its\ndescrition are sound.\nOverall, and I recommend publication in Nature Communications after the following points have been adequately\naddressed:\n- It would be helpful in Fig. 1b (and Fig. 3a) to indicate the Fermi level. It’s hard to follow what the THz field does, and what is the “dark” energy level alignment (without THz) and the one at VTHz,pk. Also, I assume the defect levels align different whether the defect is charged or not? What are the level alignments relevant for charge transport from the tip into the neutral defect (Gamma_T), from the charged defect back to the tip (thin black arrow to the right), and from the charged defect into the substrate (Gamma_S)? How to those change with THz bias?\n- Page 5: “To maximize the time resolution …. the majority of LW-STM/STS experiments have so far been performed at or\nclose to Vdc = 0”: I think that’s not generally true. What is relevant is rather that only the peak of the THz half-cycle “dips”\ninto the onset of an LDOS peak to maximize time resolution. For a highly n-doped sample, this might be achieved at a bias far\nfrom Vdc = 0.\n- Fig 2c) for 1ML: it’s hard to see a second “rising edge” at V_THz,pk ~ 1.0 V (associated to LUMO+1). Why is this not more\npronounced? In the truly quasi-static LW-STM picture, I would think the LUMO+1 edge should be of similar magnitude as the\nfirst LUMO onset? But in LW-STM, at V_THz,pk = 1.2 V, the LUMO becomes occupied first by tunneling (during the rising\nedge of the half-cycle of the transient THz field), and this charge stays for a while, so I assume this affects subsequent\ntunneling into LUMO+1 within the few 100 fs wide THz half-cycle, and hence the edge of LUMO+1 tunneling?\n- How do the (dQ_LW)/(dV_THz) curves of Fig2c look like? Do they help revealing the onsets more clearly?\n- Is Figure 2d (top left) a LW-STM map of the LUMO only, or rather the LUMO+1 (even if the orbitals look the same)? What\nare the relative contributions? Do the LW-STM images look different at V_THz,pk = 0.7 V and at V_THz,pk = 1.2 V?\n- Why does the orbital in Fig. 3b) on 2ML look different than the orbital in STS and LW-STM in Fig. 2?\n- For larger V_DC close to V_LUMO, back tunneling is effectively suppressed due to FC blockade and orbital contrast is “re-\nestablished”. How do the Q_LW-V_pk curves look under these conditions? In particular, do they reveal onsets of orbital\ntunneling? In principle this is the data shown in FigS4 e,d. There are no onsets visible, why not?\n- Plotting Q_LW versus the total bias in Fig. S4 is confusing, as without carefully checking the text it is not clear what is\nvaried (Vdc or VTHz). I recommend plotting versus VTHz.\n- The authors write “the strong suppression of back tunneling at elevated bias is attributed to a reduction in the Franck–\nCondon matrix elements.” Franck–Condon factors are strictly geometrical (dependent on displacement and phonon coupling\nstrength). What exactly causes them to vary with Vdc? Isn’t it more accurate to say that the number of energetically allowed\ntransitions changes with Vdc, rather than the Franck–Condon factors themselves\n- What determines the orbital contrast in the images Fig. 4 as a function of tip–sample distance? To what extent is the\nobserved spatial signal a convolution of the defect orbital (LDOS) and the local variation in Franck–Condon blockade\nefficiency or tunneling asymmetry? Could local variations in tip coupling or vibrational relaxation affect the apparent shape?\nThe authors explain later in the context of Fig. 5 that ΔV depends on location, but it would help to discuss this already earlier\nin the context of Fig. 4.\n- How robust are the conclusions about suppression of back tunneling drawn from the model against the greatly simplified\nassumption of a single vibrational mode?\n- In a recent publication by Roelcke et al, Nat. Photonics (2024), it is shown that THz-pulses can induce coherent lattice\nvibrations on WSe2/Au heterostructures, modulating defect levels by few meV. Why are such modulations absent in the\nWSe2/graphene systems? Are they sometimes observed? I assume this could potentially modulate the FC-blockade as well?\n- Prior studies indicate that LUMO and LUMO+1 states exhibit different charge-state lifetimes. Given the high sensitivity of\nτ_eff to changes in ΔV and Z, could these lifetimes be discerned by choosing appropriate (ΔV,Z_0), such that either state\ndominates the response? For example, since the Q_LW - Z traces can be related to different tunneling rates in the first\nbarrier, have you observed cases with a double maxima feature that might correspond to the LUMO+1 and LUMO dominated\nregimes? Overall, the role of LUMO+1 could be more discussed in the manuscript.\n- Could local variations in the graphene Fermi level (due to residual doping or charge puddles) shift the defect level\nalignment, suppressing tunneling at certain ΔV? How do local potential variations affects the observed FC blockade and\nrectification behavior at the atomic scale?\n- Are the authors able to discern between the top layer and bottom layer defect configuration lifetimes in the 2ML? I\nunderstand that the difference is small from saturation measurements but isn’t it surprising giving the different distance to the\ngraphene?\nMinor comments:\n- In the sentence where references 21,22 are cited for the first time, I am surprised to see an trARPES review, as this does\nnot provide pm spatial resolution. Do the authors refer to (ultrafast) orbital tomography using momentum microscopy (which\nis not discussed in this review)?\n- The last sentence of the second paragraph should say “…observation of ultrafast single electron dynamics…”, making\nclear the investigated regime is significantly different from previous work looking for example at hot carriers photoexcited into\ndelocalized continuum states.\n- The idea of adding an equivalent circuit diagram in Figure 1c) is nice, but it is not really clear how that helps the discussion.\nWhat is the additional information retrieved from it? This is not explained. I recommend to use this space for more detailed\njunction diagrams.\n- Page 4: in this discussion, I think it helps to refer to “orbitals” rather than “states” to clearly distinguish the scenario of\ntunneling into continuum or many degenerate states, where many electrons can tunnel at once.\n- A short note that the schemes in Figures 3d) and 3e) discuss the situation of the defect being in the charged state (so after\nTHz excitation) would be helpful (as not excitation is indicated in the scheme. Maybe add an “excitation”?).", "reviewer_3": "In their work, the authors observe ultrafast Coulomb blockade in atomic defects, specifically selenium vacancies in WSe₂\nmonolayers and bilayers. They demonstrate the capability to resolve population dynamics on the picosecond timescale and\npicometer spatial scale using lightwave-driven STM. Moreover, they achieve precise control over the states induced by\nsingle-cycle THz pulses resonant with desired transitions between the defect (vacancy) states.\nOverall, this work is of broad interest as it combines ultrafast measurement techniques with complex effects, such as Franck-\nCondon blockade, to study charge dynamics in single-atom quantum dots within TMD monolayers and bilayers. The authors\nalso investigate the optimal tip-to-defect distance, balancing enhanced coupling during the short duty cycle with the need to\nsuppress back tunneling, which otherwise reduces the probability of an electron reaching the drain\n(substrate). Franck-\nCondon blockade, a key mechanism in this study, suppresses the relaxation of the charged defect state back to neutral by\ninjecting an electron back into the STM tip. The mechanism is effective when the tip chemical potential is sufficiently high\n(close to the lowest unoccupied orbital of the defect). This creates an effective valve for charge transport, allowing the\nauthors to achieve a regime of strong STM-defect coupling during the control pulse and relatively slow defect relaxation via\nelectron emission to the drain. This setup provides high temporal and spacial resolution for monitoring defect’s dynamics\nusing delayed probe pulses. Notably, in the case of bilayers (2ML), the relaxation time is an order of magnitude longer than\nin the monolayer case, suggesting a fundamentally different transport regime between the defect and the substrate. The\nauthors propose that this may be due to suppressed tunneling through the second WSe₂ layer. These conclusions are\nsupported by a simple theoretical model that captures the tunneling processes across different coupling regimes. To\ndescribe Franck-Condon effect and its dependence on the DC bias, the authors applied a model consisting of a harmonic\npotential with a single vibrational mode for both the charged and neutral states of the defect. Then they describe the\ncharging and discharging of the defect with a Master equation. For the simulation the authors used the parameters extracted\nfrom the static measurements of the selenium vacancies. Overall, the theoretical description is adequate, taking into account\nall the states that can play a significant role in the observed process in a coherent way.\nIn addition to the general considerations, I have a more specific question for the authors. When you study the 2ML WSe₂\nsystem, is this a 2H-stacked bilayer? This configuration is naturally occurring and was extensively investigated in early\nstudies of WSe₂ bilayers. However, more recent experiments have explored different relative orientations of the layers,\nwhich can significantly alter their properties. For example, some stacking configurations do not strongly suppress single-\nelectron tunneling between the layers, as they allow for tunneling between valleys with the same spin splitting sign. This\ncould, in turn, influence the electron tunneling from selenium vacancies to the drain through the second WSe₂ layer.\nTherefore, I would like to ask if you have information on the stacking configuration in your setup. It would be interesting to\nunderstand the exact tunneling mechanism in the 2ML case, whether it involves direct tunneling through the second layer or\na more complex virtual process involving both layers. However, I understand that this question may be beyond the current\nscope of the manuscript.\nThe manuscript is well-written and comprehensive. The authors clearly introduce the basic concept of charge dynamics\nbetween the lightwave-driven STM, defect, and drain, followed by a detailed discussion of the technical aspects and distinct\nregimes in both single-layer and bilayer systems. The results are of both fundamental interest and present an important\ntechnical advancement in ultrafast measurement techniques, notably by introducing a backflow valve through the Franck-\nCondon blockade mechanism. Therefore, I recommend the manuscript for publication in Nature Communications.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts."}, "author_response": "Response to Referee Report: NCOMMS-25-31428-\nT\nWe thank all referees for their very positive and constructive feedback to our work. In the\nfollowing, we provide a point-by-point reply to all comments and questions.\nReviewer: 1\nReviewer's comment. Allerbeck et al. report the lightwave-driven manipulation of a charge state\nof a single atomic vacancy in mono- and bilayer WSe2. Building on a Franck-Condon suppression\nof back-tunnelling and lightwave-driven STM, they develop a new method that allows tracking\nthe occupation dynamics on atomic scales and with ps to fs temporal resolution. By probing the\nvacancy faster than a characteristic lifetime of a charged state, they directly observe the Coulomb\nblockade unfolding in time. The authors build an effective dynamical model to describe the role\nof tip-vacancy coupling strength in the Coulomb blockade measurements. The experiments are\nperformed at the very high technical level, and the data supports the main claims of the\nmanuscript.\nBy merging the atomic-femtosecond resolution of the lightwave-driven microscopy with the\ndirect observation of Coulomb blockade, this manuscript provides pivotal insights to the long-\nstanding issue of directly measuring the charge state dynamics in single molecules and atomic-\nscale building blocks of modern nanotechnology. Mapping the dynamic effects of coupling\nstrength provides unprecedented information for designing lightwave electronics. This\nmanuscript has the potential to become a journal highlight, and I am happy to recommend it for\npublication in Nature Communications after a minor revision.\nAuthors' answer. We thank the reviewer for their positive and thoughtful assessment of\nour manuscript, and for recognizing the conceptual and technical advances of our work,\nhighlighting \"pivotal insights\" to a long-standing challenge in the field.\nReviewer's comment. I would ask the authors to address the following issues:\n1. One of the key novelties of the manuscript is undoubtedly the measurements of the single\ndefect charging suppression with an ultrafast temporal resolution demonstrated in Fig. 5.\nI find the data very convincing, yet I have minor reservations about the interpretation:\n1.1 For the monolayer case, I am not yet convinced that the fit in Fig. 5b (bottom panel)\nstarting as early as 1.5 ps is justified. With DC offset bias almost aligned with the LUMO,\nit is the fine structure of the THz transients with all their trailing oscillations that\nbecomes important. Indeed, such structure is visible at least up to 2 ps in Fig. S2f and\ncan be resolved even beyond as oscillations exceeding the apparent noise level in Fig.\nS2g. I believe that (i) full cross-correlation measurements between pump and probe\npulses at the ML vacancy orbital could be shown and (ii) the effect of the trailing\noscillations quantified.\nAuthors' answer. We agree that THz waveform artifacts remain significant at early\ntime delays, particularly below 2 ps. For the monolayer (1ML) case, the Coulomb\nblockade signal itself extends over a very short interval (~2 ps), making it infeasible to\nexclude data below 2 ps without discarding essential dynamics. Although residual THz\noscillations are present at 1.5 ps (as evident in Fig. S2f,g), their influence is mitigated\nby the use of single-exponential fits, which statistically average over such\nperturbations. This procedure results in increased uncertainty rather than systematic\nbias.\nAction taken: We now explicitly state in the discussion of Fig. 5: \"We determine τeff\nvia exponential fitting of the spectra at time delays between 1.5 ps and 12 ps, ensuring\nrobust convergence of fit parameters. The error bar includes disturbances related to\nTHz trailing oscillations at high V .\"\ndc\n1.2 Reviewer's comment. What is the effect of stochastic tunnelling on the ultrafast\nCoulomb blockade measurements when V approaches V (e.g. in Fig. 5c)? In Fig.\ndc LUMO\nS1, the authors clearly show that at z = z – 3.5 Å and V ~ V + 0.2V the DC current\n0 dc LUMO\nshows Coulomb blockade signatures in a purely static experiment. How does the\nstochastic occupation compare to the lightwave-driven one? I believe that Fig. 5c\nwould benefit from showing a DC current curve as a function of V detuning.\ndc\nAuthors' answer. In all dynamic measurements, V is kept intentionally below the\ndc\nthreshold for stochastic tunneling. Already at V - V ~ 0.02V, quasi-static\ndc LUMO\noccupation from stochastic tunneling events dominates over any dynamic signature,\nwhich results in increased noise levels in the dynamic measurements, as shown in the\nbottom panel of Fig. 3c and third panel of Fig. S5c. The breakdown of lightwave-driven\nsignals due to a quasi-static population of stochastic tunneling events shown in Fig. 3c\nis well reproduced by our model simulation. Hence, transient measurements at V ~\ndc\nV + 0.2V are not feasible.\nLUMO\nBelow we present the dc current as we detune the voltage for the dynamic\nmeasurements in Fig. 5c. The standard deviation on the dc current during the dynamic\nmeasurements exceeds the magnitude of the dynamic signal already for a detuning\n<-10 mV.\nAction taken: In the discussion of Fig.5 c,f we added: \"The onset of quasi-static\ncharging at dc biases above the LUMO resonance, as shown in Fig. 3c, precludes\nultrafast measurements due to dominant stochastic tunneling.\"\n1.3 Reviewer's comment. Semi-log scale in exponential decay curves would be more\nconvincing and would allow the reader better estimate the noise level and fit reliability.\nAuthors' answer. A semi-logarithmic representation would indeed be helpful to\nvisualize data that covers multiple orders of magnitude. In our case, we observe a\nsuppression of 50% which is – in our opinion – optimally presented in a linear scale.\nMoreover, the correct presentation of THz reflections and artifacts from waveform\nsampling (linear effect) benefits from a linear representation. Offsets used in the data\nrepresentation, such as Q eq, further complicate a correct representation on a semi-\nLW\nlog scale. For these reasons, we suggest to keep the presentation in the current form.\n2. Reviewer's comment. Is it possible to image the vacancy orbital in 2ML with LW-STM\n(e.g. in Fig. 2d (iii)) by increasing the integration time? In Fig. 2c, with a larger integration\ntime, the turquoise curve is above the grey one, suggesting that the orbital should be\nvisible as compared to the surroundings. If yes (e.g. a linecut across with orbital with 100x\nintegration time), could this estimate be used to quantify the back-tunnelling rate?\nAuthors' answer. We acknowledge that the slight spectral difference in Fig. 2c suggests\na marginal orbital signal even for Vac /2ML using V = 0 V, where the contrast due to\nSe dc\nback tunneling is extremely low. However, in practice, increased integration times did not\nresult in discernible orbital contrast in imaging conditions. Due to technical limitations\nsuch as drift, the integration time cannot be easily increased by a factor of 100x since these\nLW-STM images must be taken without feedback.\n3. Reviewer's comment. The authors discuss possible vacancy level shifts due to coherent\nlattice dynamics introduced by the THz pump pulse, comparing their work to Ref. 28. Yet,\nthe present data shows no out-of-plane coherent lattice motion in either the pristine ML\nor BL, in contrast to the drum mode in Ref. 28. What kind of lattice dynamics do the authors\nimply then? In-plane? Can the authors rule out the effects of coherent acoustic phonon\nwave packets observed by Sheng et al. in PRL 129, 043001 (2022)?\nI would also take the liberty to notice, that Ref. 28 reveals dynamic shifts of the vacancy\nstates even in the absence of charging, in contrast to what is written in the present\nmanuscript.\nAuthors' answer. We appreciate the opportunity to clarify this point. The epitaxial\nsamples used in our measurements have no underlying moiré structure that promotes an\nout-of-plane coherent motion as reported in Ref. 28, hence we do not expect the same\nbehavior. Additionally, the drum-mode as observed in Ref. 28 depends on the work\nfunction of the substrate, which is significantly altered from Au(111) to graphene.\nNonetheless, inelastic electron transfer inevitably excites lattice modes. We consider both\nin-plane and weak out-of-plane phonon contributions, but fast lateral dissipation likely\nprevents their local detection. In our measurements at high V such signatures would also\ndc\nbe mixed with coherent signal modulations by multiple excitations and THz waveform\nsampling, preventing a clean distinction of both effects.\nRegarding Sheng et al., we find no evidence of substrate-related acoustic wave effects in\nour system. Presumably the inter-layer Gr/WSe heterostructure coupling is too soft for\n2\nany vertical acoustic excitation. Moreover, THz propagation inside the conductive SiC\nsubstrate is strongly attenuated, such that THz reflections remain insignificant. However,\nTHz pulse reflections from the metallic sample plate were observed when using semi-\ninsulating SiC substrates instead of conductive SiC (used here).\nAction taken: We added a paragraph discussing coherent vibrations and the findings in\nRef. 28.\n4. Reviewer's comment. I encourage the authors to compare their technique with other\nrecent works in the field of lightwave-driven microscopy. In particular, I would like to bring\nthe authors’ attention to the recent paper by Siday et al., Nature 629, 329-334 (2024) that\nreports the subcycle tunnelling currents within the THz pulse itself, and is thus sensitive to\nthe non-rectified, including the back-tunnelling, charges (e.g. Fig. 5 in the Siday et al.).\nAuthors' answer. We are aware of the potential provided by subcycle time resolution as\npresented by Siday et al. and thank the reviewer for highlighting the importance in this\ncontext. A thorough description of non-equilibrium tunneling currents within the THz\npulse is crucial for understanding the limits of net charge rectification. Our simulations\ninclude the complete THz waveform and hence accounts for the statistics of forward and\nbackward tunneling along the transient profile of the THz field. In contrast to the work of\nSiday et al., our lightwave-driven currents are limited to usually <1e/pulse and the back\ntunneling process is in our case not coherently driven by the pulse but underlies the\nstochastic tunneling process.\nAction taken: We now cite the work by Sidday et al. in the conclusion paragraph.\n5. Reviewer's comment. The authors very convincingly highlight the importance of, and the\nhopes related to the defect-based functional properties, so it would be great to see what\nthe field can learn from this manuscript. Is the regime of strong coupling to the leads\naccessible in e.g. vertical heterostructures, THz or optical photonic resonators? The spatial\nprobing of the ultrafast Coulomb blockade (Figs. 4, 5) is very impressive: could it provide\nus with some ideas for the device design or be used in device simulations? Could this\ntechnique be applied for molecules? I would appreciate a broader outlook.\nAuthors' answer. Comparable coupling strengths are indeed achievable also in a device\ngeometry. Van der Waals gaps, such as used here for the bottom (sample) electrode, offers\nwell controlled coupling by the charge density (that maybe controlled by external gates)\nand spacing to adjacent layers. In a vertical tunneling device the STM tip would simply be\nreplaced by a flat electrode, separated by an ultrathin insulating spacer, such as hBN. Such\ntunneling contacts are already employed in high quality transport devices [e.g. C. Li Nature\nMaterials 22, 570 (2023)]. Facilitating tunneling channels across localized states, such as\nfrom single vacancies acting as a single electron transistor is of course much more\nchallenging using lithographically defined contacts, but potentially feasible.\nIndeed, we see many similarities between discrete in-gap states in semiconductors and\nmolecules. We anticipate that ultrafast Coulomb blockade could be probed in molecules\nadsorbed on thin insulating layers such as NaCl where the effect of the Franck-Condon\nblockade might even be stronger due to the stronger electron-phonon coupling on ionic\nsubstrates, but the required small tip-sample distances might be challenging to access for\nsome adsorbates.\n6. Reviewer's comment. Manuscript presentation suggestions:\nIs the use of the picocavity concept really relevant? The authors mention it only once and\nbriefly, in the abstract, without much context.\nAuthors' answer. In the context of this manuscript, the STM picocavity provides near-field\nenhancement of the transient voltage and spatial resolution down to the atomic level,\nhence forming it an essential ingredient for LW-STM. We believe picocavity is a commonly\nused terminology in the STM community and fits the scope of our article well. However,\nsince this terminology is not strictly required, we replaced \"picocavity\" with \"STM junction\"\nfor clarity.\nAction taken: Replace \"focused onto the picocavity of a scanning tunneling microscope\"\nto \"focused onto the junction of a scanning tunneling microscope\"\n7. Reviewer's comment. In Fig. 4, the authors seem to use the term LW-STS to describe the\nz-dependence of the signal. While it z-spectroscopy is also a kind of spectroscopy, it was\nconfusing for me as a reader.\nAuthors' answer. We understand the confusion and clarify the terminology in the caption.\nAction taken: We rephrased \"LW-STS\" in the figure caption with \"z-spectroscopy\".\nReviewer: 2\nReviewer's comment. The work by Allerbeck et al. demonstrates picosecond-resolved control\nand detection of single-electron charge dynamics in an atomic-scale quantum dot (VacSe defects\nin WSe₂) using lightwave-driven STM (LW-STM) and a THz pump–THz probe configuration. The\nkey finding is the direct observation of ultrafast Coulomb blockade and its modulation via the\nFranck–Condon (FC) blockade. I think this is a significant achievement in the field of atomic\nquantum transport.\nEven though charge state lifetimes have been studied extensively using conventional STM (see\ncited literature), the scientific results enter a new regime of real-time observation of picosecond\ncharge state lifetimes. The main novelty in this work is the ability to observe the transient Coulomb\nblockade dynamics in real time and down to a few picoseconds. In addition, the authors\ndemonstrate active control of tunneling directionality by inducing an asymmetry in forward-\nbackward tunneling by the Franck–Condon Blockade, which is a new achievement, despite the fact\nthat the underlying physics is well-established. The paper is well written, the data and figures are\npresented in good quality and the methodology and its description are sound.\nOverall, and I recommend publication in Nature Communications after the following points have\nbeen adequately addressed:\nAuthors' answer. We thank the reviewer for the constructive and thoughtful assessment.\n1. Reviewer's comment. It would be helpful in Fig. 1b (and Fig. 3a) to indicate the Fermi\nlevel. It’s hard to follow what the THz field does, and what is the “dark” energy level\nalignment (without THz) and the one at VTHz,pk. Also, I assume the defect levels align\ndifferent whether the defect is charged or not? What are the level alignments relevant for\ncharge transport from the tip into the neutral defect (Gamma_T), from the charged defect\nback to the tip (thin black arrow to the right), and from the charged defect into the\nsubstrate (Gamma_S)? How do those change with THz bias?\nAuthors' answer. The Fermi level of the sample is pinned by the QFEG substrate, as\nindicated by the filled dark gray box labeled \"QFEG\" in Fig. 1b, representing the energy\nalignment in the absence of bias and THz field (\"dark\"). The Fermi level of the substrate\nsets the Vac charge state. We could show in a previous study [Xiang et al, Nat. Comm.\nSe\n(2023)] that the defect charge state can be controlled via the substrate work function. In\nMoS , for instance, the S vacancy is negatively charged on QFEG/SiC. This Fermi level\n2\npinning is specific to the TMD/graphene heterostructure and would be different for free-\nstanding or bulk WSe . For the tip electrode the electrochemical potential can be shifted\n2\nby both a static dc bias and transient THz bias, which superimpose, indicated in Fig. 1b by\nthe light gray box and blue transient.\nThe referee raises an interesting question about the relevant defect state energies for\ntunneling in (Γ in) and out (Γ and Γ out) of the defects state. From STS measurements we\nT S T\ncan measure the electron attachment energy (charge transition energy Vac 0 → Vac -1) to\nSe Se\nthe LUMO and LUMO+1 (around 0.7 V and 0.8 V for 1ML WSe , see Fig.2a). Once the\n2\nelectron is attached to the orbital (at the ZPL or a vibrationally excited state) it will relax to\nthe vibronic ground state (cf. schematic in Fig. 3d). This relaxation energy can be anywhere\nfrom 0 to a few tens of meV resulting from a finite electron-phonon coupling visible as\nvibronic sidebands in the STS spectra. Despite the energy relaxation of the electron in and\nout of the defect state is small, it is significant enough to generate a difference in the\ntunneling rate out of the defect in the case of the tip electrode. This is the Franck-Condon\nblockade.\nThe defect states only minimally shift by a linear scaling factor with applied bias voltage\n(both dc and THz) due to the well-known double barrier tunneling junction effect.\nWhile the LUMO energy remains almost unaffected, the presence of an extra electron\nrenormalizes all other electronic states, especially the frontier orbitals. Therefore, the Vac\nSe\nstates of the charge neutral and negatively charged vacancy are indeed quite different but\ndifficult for us to assess because we cannot independently set the charge state and\nperform STS. Importantly for the current manuscript, the electron attachment energy for\nthe transition from Vac -1 → Vac -2 is shifted outside the band gap because of strong\nSe Se\nCoulomb repulsion, resulting in the observed Coulomb blockade.\nAction taken: We added a label for the sample Fermi level in Fig. 1b. We mention that the\nFermi level is pinned by the QFEG substrate in the caption of Fig. 1 and the Methods\nsection.\n2. Reviewer's comment. Page 5: “To maximize the time resolution …. the majority of LW-\nSTM/STS experiments have so far been performed at or close to Vdc = 0”: I think that’s\nnot generally true. What is relevant is rather that only the peak of the THz half-cycle “dips”\ninto the onset of an LDOS peak to maximize time resolution. For a highly n-doped sample,\nthis might be achieved at a bias far from Vdc = 0.\nAuthors' answer. The reviewer correctly notes that time resolution is improved when only\nthe very peak of the THz field interacts with the orbital state, which may be enhanced by\nshifting V far from the resonance. Relevant studies [Refs. 23,27,28,41] predominantly\ndc\nused V ≪ V pk. On metallic substrates a small V facilitates active feedback, however,\ndc THz dc\nstatic currents from dc tunneling (background) weaken the sensitivity of dynamics\nmeasurements for higher dc voltage.\nAction taken: We clarified the sentence on time resolution.\n3. Reviewer's comment. Fig 2c) for 1ML: it’s hard to see a second “rising edge” at V_THz,pk\n~ 1.0 V (associated to LUMO+1). Why is this not more pronounced? In the truly quasi-\nstatic LW-STM picture, I would think the LUMO+1 edge should be of similar magnitude as\nthe first LUMO onset? But in LW-STM, at V_THz,pk = 1.2 V, the LUMO becomes occupied\nfirst by tunneling (during the rising edge of the half-cycle of the transient THz field), and\nthis charge stays for a while, so I assume this affects subsequent tunneling into LUMO+1\nwithin the few 100 fs wide THz half-cycle, and hence the edge of LUMO+1 tunneling?\nAuthors' answer. We thank the reviewer for this question. What the reviewer notes as an\napparent reduced amplitude of the LUMO+1 onset is mainly a consequence of the\nexponential scaling of the plot (discussed below). In the numerical derivative (dQ )/(dV )\nLW THz\n(see next question) both states yield a comparable increase in the LW tunneling current.\nHowever, indeed V pk voltage sweeps and V sweeps are not equivalent since increasing\nTHz dc\nV pk will always slightly increase the pulse duration (duty cycle) of the lower lying states,\nTHz\nhence gradually increasing their contribution to the rectified charge. Accordingly, we\nobserve a slope in Fig. 2c and not a step like increase with each state as in I(V)\nmeasurements.\nThere is an additional effect related to current density: If voltage sweeps are performed\nclose to current saturation conditions (e.g. Q ~1 e/pulse for charge state lifetimes longer\nLW\nthan the pulse duration or equivalently I~I in dc STS as reported in Ref. 16) then the\nsat\nLUMO/LUMO+1 tunneling channels do not contribute independently and start to block\neach other. Therefore, when hitting the LUMO+1 threshold the additional I(V)/ Q (V )\nLW THz\nstep might be reduced or fully quenched.\nThese factors can explain the differences between Q(V_THz) and I(V) curves. The quasi-\nstatic LW-STM picture still applies.\n4. Reviewer's comment. How do the (dQ_LW)/(dV_THz) curves of Fig. 2c look like? Do they\nhelp revealing the onsets more clearly?\nAuthors' answer.\nThe figure below shows the corresponding (dQ )/(dV ) curve. Indeed, LUMO and\nLW THz\nLUMO+1 show similar amplitude, as expected in a regime of Q ≪1 e/pulse.\nLW\nHowever,\nowing to low sampling resolution of Q - V pk curves in Fig. 2c, which is beneficial for\nLW THz\nmeasurement stability, the quality of the numeric derivative is low, and we decided that\nthere is no significant benefit of showing this additional plot in the manuscript.\n.\n5. Reviewer's comment. Is Figure 2d (top left) a LW-STM map of the LUMO only, or rather\nthe LUMO+1 (even if the orbitals look the same)? What are the relative contributions? Do\nthe LW-STM images look different at V_THz,pk = 0.7 V and at V_THz,pk = 1.2 V?\nAuthors' answer. The orbital in Fig. 2d(i) is predominantly the LUMO with contributions\nof the LUMO+1 for the same arguments as discussed above. The orbital shape of LUMO\nand LUMO+1 is very similar (Fig. S1) and indistinguishable with LW-STM contrast of\npresent experiments. Also, in dc-STM measurements the LUMO and LUMO+1 are not\nsignificantly different, cf. Ref 18.\n6. Reviewer's comment. Why does the orbital in Fig. 3b) on 2ML look different than the\norbital in STS and LW-STM in Fig. 2?\nAuthors' answer. We thank the referee for this keen observation. The orbital in Fig. 3b is\nrecorded at high V , which is required to reduce back tunneling on 2ML. At a dc bias close\ndc\nto orbital, the LW-STM orbital shape can become saturated due to effective dynamic\ncharging of the orbital as observed with dc STM in Ref. 16 Fig. 3c. In contrast, Vac /1ML\nSe\nin Fig. 2d(i) acquired at low V appears unsaturated due to fast charge relaxation much\ndc\nshorter than the pulse duration. Fig. 4 shows the transition from an unsaturated to the\nsaturated orbital of Vac /2ML. There is also a height dependence on the shape of the\nSe\norbital as shown in ref 16. In Fig. 3b we observe a very similar orbital shape for Vac /2ML\nSe\nfor large tip-sample distances as in the 1ML case.\nAction taken: For clarity, we added a comment about the orbital shape in the discussion\nof Fig. 3b. \"Note that the orbital shape of Vac in Fig. 3b appears saturated, unlike in the\nSe\nmonolayer case in Fig. 2d(i). Due to the required ΔV ~ 0V, dc electrons block the THz\nelectrons at some tip positions. The same orbital shape as for the 1ML case can be\nrecovered for specific V and tip heights as shown in Fig. 4b\"\ndc\n7. Reviewer's comment. For larger V_DC close to V_LUMO, back tunneling is effectively\nsuppressed due to FC blockade and orbital contrast is “re-established”. How do the Q_LW-\nV_pk curves look under these conditions? In particular, do they reveal onsets of orbital\ntunneling? In principle this is the data shown in FigS4 e,d. There are no onsets visible, why\nnot?\nAuthors' answer. We refer to Fig. S4 b,e where Q -V pk curves for different V are\nLW THz dc\nshown. For the bilayer in Fig. S4e, by approaching the LUMO resonance with the dc bias:\n700 mV (dark blue) --> 730 mV (green), the charge rectification is gradually increased.\nIrrespective of the dc bias, orbital tunneling sets in at the same V + V pk value\ndc THz\ncorresponding to resonant tunneling into the LUMO as shown by the black curve in the\nsame plot.\nWhile the FC-blockade effectively suppresses back tunneling from LUMO, V remains far\ndc\nfrom the onset of LUMO+1, hence the signal is dominated by LUMO and no additional\nonset for the LUMO+1 orbital is observed.\nAction taken: We added a comment in the caption of Fig. S4 for better understanding.\n8. Reviewer's comment. Plotting Q versus the total bias in Fig. S4 is confusing, as without\nLW\ncarefully checking the text it is not clear what varies (V or V ). I recommend plotting\ndc THz\nversus V .\nTHz\nAuthors' answer. We understand that the axis label may be unintuitive without carefully\nreading the figure caption.\nAction taken: We improved the axis label in the figure to V pk (V) offset by V .\nTHz dc\n9. Reviewer's comment. The authors write “the strong suppression of back tunneling at\nelevated bias is attributed to a reduction in the Franck–Condon matrix elements.” Franck–\nCondon factors are strictly geometrical (dependent on displacement and phonon coupling\nstrength). What exactly causes them to vary with Vdc? Isn’t it more accurate to say that the\nnumber of energetically allowed transitions changes with Vdc, rather than the Franck–\nCondon factors themselves?\nAuthors' answer. We thank the referee for pointing out this mistake. The Franck-Condon\nfactors are in good approximation bias independent, assuming the displacement does not\nvary significantly with the applied electric field in the junction. The suppression of the back\ntunneling described here is rather caused by the reduction in the number of energetically\nallowed transitions, as the referee pointed out correctly.\nAction taken: We rewrote the sentence to \"Importantly, the strong suppression of back\ntunneling at elevated bias is attributed to a reduction in the energetically allowed vibronic\ntransitions into the neutral state.\".\n10. Reviewer's comment. What determines the orbital contrast in the images Fig. 4 as a\nfunction of tip–sample distance? To what extent is the observed spatial signal a\nconvolution of the defect orbital (LDOS) and the local variation in Franck–Condon\nblockade efficiency or tunneling asymmetry? Could local variations in tip coupling or\nvibrational relaxation affect the apparent shape? The authors explain later in the context\nof Fig. 5 that ΔV depends on location, but it would help to discuss this already earlier in\nthe context of Fig. 4.\nAuthors' answer. The reviewer correctly notes that the images show a convolution of\nLDOS, Coulomb and Franck-Condon blockade. The series of orbital images in Fig. 4\nrepresents the transition from an unsaturated shape (dz = -1 Å) to the saturated orbital\n(dz = -3 Å), related to an increased Coulomb blocking effect from dc electrons. At dz = -\n3.5 Å, charge state depletion by the tip (back tunneling) dominates over tunneling\nchannels to the substrate resulting in net zero charge rectification close to the defect\ncenter.\nAction taken: We added a sentence to the discussion of Fig. 3b.\n11. Reviewer's comment. How robust are the conclusions about suppression of back\ntunneling drawn from the model against the greatly simplified assumption of a single\nvibrational mode?\nAuthors' answer. The simplification with a single vibrational mode is robust and\nsufficiently captures the relevant processes for bias voltages several meV below the LUMO\nonset. The simplified model breaks down as the bias voltage reaches E since it does\nLUMO\nnot capture the additional Franck-Condon blockade due to the low-energy vibrational\nmodes. However, the general conclusions about the suppression of back tunneling is still\nrobust, as the effect of the low-energy modes becomes negligible for sufficiently large\nenergy differences. Adding a second vibrational mode will lead to unreasonable\ncomputational costs with minor accuracy gain for most energy ranges. We agree that\nfurther vibrational modes can account for minor differences between simulation and\nexperiment as V approaches E .\ndc LUMO\n12. Reviewer's comment. In a recent publication by Roelcke et al, Nat. Photonics (2024), it is\nshown that THz-pulses can induce coherent lattice vibrations on WSe /Au\n2\nheterostructures, modulating defect levels by few meV. Why are such modulations absent\nin the WSe /graphene systems? Are they sometimes observed? I assume this could\n2\npotentially modulate the FC-blockade as well?\nAuthors' answer. We would expect that an energy modulation of the defect state within\na few 10 mV driven by a coherent motion of the WSe layer as investigated by Roelcke et\n2\nal. to also modulate the Franck–Condon blockade significantly. The lack of a WSe\n2\nmoiré\nstructure in our sample reduces mechanical strain and may facilitate fast delocalization of\nvibrational excitations. At elevated V\ndc\nthe decoherence of lattice excitations is amplified\nby charge transfers from pulse replicas. Neither on bare WSe nor on the Vac orbitals, we\nobserve a modulation related to vibrational excitation or energetic shifts of the defect\nlevels that are significant over the modulation by the THz pulse train.\nAction taken: We added a discussion on coherent signal modulation in our experiments.\n13. Reviewer's comment. Prior studies indicate that LUMO and LUMO+1 states exhibit\ndifferent charge-state lifetimes. Given the high sensitivity of τ to changes in ΔV and ,\neff z\ncould these lifetimes be discerned by choosing appropriate (ΔV,z ), such that either state\n0\ndominates the response? For example, since the Q – z traces can be related to different\nLW\ntunneling rates in the first barrier, have you observed cases with a double maxima feature\nthat might correspond to the LUMO+1 and LUMO dominated regimes? Overall, the role\nof LUMO+1 could be more discussed in the manuscript.\nAuthors' answer. In our transient experiments ΔV is much smaller than the energy\nseparation of the LUMO and LUMO+1, hence the FC-blockade is ineffective for the\nLUMO+1. At the same time, V must remain below the LUMO resonance to avoid large\ndc\ndc currents. For this reason, we cannot effectively measure the LUMO+1 independently.\nPrior studies (Ref. 16) were able to extract the charge-state lifetime of the LUMO+1 only\nfrom an effective measurement of both defect levels as well as of the LUMO alone.\n14. Reviewer's comment. Could local variations in the graphene Fermi level (due to residual\ndoping or charge puddles) shift the defect level alignment, suppressing tunneling at\ncertain ΔV? How do local potential variations affects the observed FC blockade and\nrectification behavior at the atomic scale?\nAuthors' answer. The Fermi level alignment in our samples is set by the QFEG substrate,\nwhere unintentionally charged defects in the WSe layer locally shift the local potential.\nSe 2\nvariation at the meV level remains. For this reason, we define V differently for every\nLUMO\nVac , such that the relative alignment ΔV is comparable within the margin of error. Local\nSe\nvariations of the potential and defect complexes, e.g. Vac +Mo , impact the charge-state\nSe W\nlifetime significantly as discussed in Ref. 16 and seen for different Vac in this study (Fig.\nSe\nS1). Fluctuations from the QFEG itself are not observed due to the high quality of the\nsubstrate on the length scale of our measurements.\n15. Reviewer's comment. Are the authors able to discern between the top layer and bottom\nlayer defect configuration lifetimes in the 2ML? I understand that the difference is small\nfrom saturation measurements but isn’t it surprising giving the different distance to the\ngraphene?\nAuthors' answer. The STM enables a clear distinction between top or bottom Vac in\nSe\nupper or lower layer of bilayer WSe , as studied in Ref. 16, Fig. S5 and Fig. S8. In agreement\n2\nwith our previous report, we confirm that the top and bottom Vac in the same WSe layer\nSe 2\nhave a very similar charge state lifetime. While structurally the separation from the\ngraphene varies between the two cases, the vertical distribution of the electronic defect\norbital within the layer is actually quite similar. For this reason, the comparable tunneling\nrates are in fact not that surprising. As expected, the charge state lifetime of Vac in\nSe\ndifferent layers is substantially different. In the present study we only investigated Vac in\nSe\nthe upper WSe layer, where lower layers provide decoupling from the substrate (Fig. S1).\n2\n16. Reviewer's comment. Minor comments:\nIn the sentence where references 21,22 are cited for the first time, I am surprised to see an\ntrARPES review, as this does not provide pm spatial resolution. Do the authors refer to\n(ultrafast) orbital tomography using momentum microscopy (which is not discussed in this\nreview)?\nAuthors' answer. We agree with the reviewer that a reference to (ultrafast) orbital\ntomography is a better fit to support this sentence.\nAction taken: We cite R. Wallauer et al. Science 371, 1056 (2021) instead.\n17. Reviewer's comment. The last sentence of the second paragraph should say\n“…observation of ultrafast single electron dynamics…”, making clear the investigated\nregime is significantly different from previous work looking for example at hot carriers\nphotoexcited into delocalized continuum states.\nAuthors' answer. We thank the reviewer for this valuable comment.\nAction taken: We adapted the sentence accordingly.\n18. Reviewer's comment. The idea of adding an equivalent circuit diagram in Figure 1c) is\nnice, but it is not really clear how that helps the discussion. What is the additional\ninformation retrieved from it? This is not explained. I recommend to use this space for\nmore detailed junction diagrams.\nAuthors' answer. We agree with the reviewer that the equivalent circuit diagram does\nnot necessarily provide additional information. Instead, it is intended as an interdisciplinary\nbridge, connecting typical visualizations in the STM and quantum dot community. This\naligns well with previous illustrations, e.g. in T. Aliyar et al. Nano Lett. 24, 2142 (2024), and\nwe insist to keep it to highlight the interdisciplinary aspect of this study.\nThe junction diagram is kept as simple as possible to address a broad audience, while\ncontaining all relevant information. Additional diagrams in Fig. 3, Fig. 4, and Fig. S6\nelaborate the details.\n19. Reviewer's comment. Page 4: in this discussion, I think it helps to refer to “orbitals” rather\nthan “states” to clearly distinguish the scenario of tunneling into continuum or many\ndegenerate states, where many electrons can tunnel at once.\nAuthors' answer. We understand that the terminology may have been misleading here.\nIn scope of the introduction, we would like address excitations and states in a general\ncontext, although the Vac localized states refer to orbitals.\nSe\nAction taken: We made clear that we refer to \"localized\" states (orbitals) at this part of\nthe discussion.\n20. Reviewer's comment. A short note that the schemes in Figures 3d) and 3e) discuss the\nsituation of the defect being in the charged state (so after THz excitation) would be helpful\n(as not excitation is indicated in the scheme. Maybe add an “excitation”?).\nAuthors' answer. The schematics in Figs. 3d and 3e show a transition from charged to\nneutral state, which in context of the article requires an excitation. To ease the\nunderstanding from the figure alone we agree to adapt the labels.\nAction taken: In the figure we adapted the labelling to “charged (excited)”.\nReviewer: 3\nReviewer's comment. In their work, the authors observe ultrafast Coulomb blockade in atomic\ndefects, specifically selenium vacancies in WSe₂ monolayers and bilayers. They demonstrate the\ncapability to resolve population dynamics on the picosecond timescale and picometer spatial\nscale using lightwave-driven STM. Moreover, they achieve precise control over the states induced\nby single-cycle THz pulses resonant with desired transitions between the defect (vacancy) states.\nOverall, this work is of broad interest as it combines ultrafast measurement techniques with\ncomplex effects, such as Franck-Condon blockade, to study charge dynamics in single-atom\nquantum dots within TMD monolayers and bilayers. The authors also investigate the optimal tip-\nto-defect distance, balancing enhanced coupling during the short duty cycle with the need to\nsuppress back tunneling, which otherwise reduces the probability of an electron reaching the drain\n(substrate). Franck-Condon blockade, a key mechanism in this study, suppresses the relaxation of\nthe charged defect state back to neutral by injecting an electron back into the STM tip. The\nmechanism is effective when the tip chemical potential is sufficiently high (close to the lowest\nunoccupied orbital of the defect). This creates an effective valve for charge transport, allowing the\nauthors to achieve a regime of strong STM-defect coupling during the control pulse and relatively\nslow defect relaxation via electron emission to the drain. This setup provides high temporal and\nspatial resolution for monitoring defect’s dynamics using delayed probe pulses. Notably, in the\ncase of bilayers (2ML), the relaxation time is an order of magnitude longer than in the monolayer\ncase, suggesting a fundamentally different transport regime between the defect and the substrate.\nThe authors propose that this may be due to suppressed tunneling through the second WSe₂\nlayer. These conclusions are supported by a simple theoretical model that captures the tunneling\nprocesses across different coupling regimes. To describe Franck-Condon effect and its\ndependence on the DC bias, the authors applied a model consisting of a harmonic potential with\na single vibrational mode for both the charged and neutral states of the defect. Then they describe\nthe charging and discharging of the defect with a Master equation. For the simulation the authors\nused the parameters extracted from the static measurements of the selenium vacancies. Overall,\nthe theoretical description is adequate, taking into account all the states that can play a significant\nrole in the observed process in a coherent way.\nAuthors' answer. We thank the reviewer for their positive and constructive feedback, calling the\nwork \"of broad interest\".\n1. Reviewer's comment. In addition to the general considerations, I have a more specific\nquestion for the authors. When you study the 2ML WSe₂ system, is this a 2H-stacked\nbilayer? This configuration is naturally occurring and was extensively investigated in early\nstudies of WSe₂ bilayers. However, more recent experiments have explored different\nrelative orientations of the layers, which can significantly alter their properties. For\nexample, some stacking configurations do not strongly suppress single-electron tunneling\nbetween the layers, as they allow for tunneling between valleys with the same spin splitting\nsign. This could, in turn, influence the electron tunneling from selenium vacancies to the\ndrain through the second WSe₂ layer. Therefore, I would like to ask if you have information\non the stacking configuration in your setup.\nAuthors' answer. We thank the referee for raising this interesting point. Our epitaxial\nWSe samples exhibit varying stacking configurations, including 2H and 3R, depending on\n2\nthe specific island. The stacking order can, in principle, be inferred from large-area STM\nimages from the relative orientation of the triangular islands, or the orientation of the\nthree-fold symmetric defect orbitals that are either equivalently oriented or mirrored.\nWhile variations in interlayer spacing due to stacking may influence the charge-state\nlifetime, we expect their effect to be secondary compared to fluctuations in the local\nelectrostatic environment, such as those arising from nearby charged defects. See also Ref.\n16 for more elaborate statistics. We did not identify the stacking order of the bilayer region\nused for the present study.\n2. Reviewer's comment. It would be interesting to understand the exact tunneling\nmechanism in the 2ML case, whether it involves direct tunneling through the second layer\nor a more complex virtual process involving both layers. However, I understand that this\nquestion may be beyond the current scope of the manuscript.\nAuthors' answer. As studied in Ref. 16 (Fig. 4) the lower layer does indeed play a role in\nthe tunneling process since the Vac defect orbitals are partially delocalized to the lower\nSe\nlayer due to hybridization. Therefore, the charge-state lifetime of Vac / WSe increases\nSe 2\nsub-exponentially with additional layers due finite orbital density of states in adjacent\nlayers, in contrast to molecules on insulating NaCl decoupling layers [see Kaiser et al,\n2020]. The layer-dependent scaling of the charge-state lifetime can be understood by first\norder tunneling processes.\nReviewer's comment. The manuscript is well-written and comprehensive. The authors clearly\nintroduce the basic concept of charge dynamics between the lightwave-driven STM, defect, and\ndrain, followed by a detailed discussion of the technical aspects and distinct regimes in both\nsingle-layer and bilayer systems. The results are of both fundamental interest and present an\nimportant technical advancement in ultrafast measurement techniques, notably by introducing a\nbackflow valve through the Franck-Condon blockade mechanism. Therefore, I recommend the\nmanuscript for publication in Nature Communications.\nAuthors' answer. We thank the reviewer once again for their thoughtful feedback and their\nrecommendation for publication.\nResponse to Referee Report: NCOMMS-25-31428A\nWe thank all referees for their very positive and constructive feedback to our work. In the\nfollowing, we provide a point-by-point reply to all comments and questions.\nReviewer: 1\nReviewer's comment. I am grateful to the authors for their insightful responses to the\nquestions raised in the review. The revised manuscript addresses all the issues, and I am happy\nto recommend it for publication in Nature Communications.\nAuthors' answer. We thank the reviewer for the scientific discussion, and we are\ngrateful for the positive assessment.\nReviewer: 2\nReviewer's comment. I thank the authors for their detailed and instructive response to the\npoints raised. I am pleased to see that all questions are clarified and that the authors addressed\nthe comments in a satisfactory manner. The revisions have improved the clarity and quality of\nthe manuscript, and I do not have any major concern.\nTherefore, I recommend the manuscript for publication in Nature Communications.\nAuthors' answer. We thank the reviewer for the scientific discussion and we are grateful\nfor the positive assessment.\nReviewer's comment. I have one minor follow-up question to one of the comments, for which I\nwould be happy to receive an answer:\n12. Reviewer's comment: In a recent publication by Roelcke et al, Nat. Photonics (2024), it is\nshown that THz-pulses can induce coherent lattice vibrations on WSe2/Au heterostructures,\nmodulating defect levels by few meV. Why are such modulations absent in the WSe2/graphene\nsystems? Are they sometimes observed? I assume this could potentially modulate the FC-\nblockade as well?\nAuthors' answer. We would expect that an energy modulation of the defect state within\na few 10 mV driven by a coherent motion of the WSe2 layer as investigated by Roelcke\net al. to also modulate the Franck–Condon blockade significantly. The lack of a WSe2\nmoiré structure in our sample reduces mechanical strain and may facilitate fast\ndelocalization of vibrational excitations. At elevated V the decoherence of lattice\ndc\nexcitations is amplified by charge transfers from pulse replicas. Neither on bare WSe2\nnor on the VacSe orbitals, we observe a modulation related to vibrational excitation or\nenergetic shifts of the defect levels that are significant over the modulation by the THz\npulse train.\nAction taken: We added a discussion on coherent signal modulation in our experiments\nReviewer's comment. Follow-up question: Thank you for the clarification. I appreciate that the\nabsence of a moiré pattern and associated mechanical strain may lead to faster delocalization of\nvibrational excitations, potentially reducing the amplitude or coherence of localized modes.\nNeverthesless, it is less clear to me why it should be a consequence that no drum mode is\nobserved at all, as the authors write in their manuscript. I would expect such a mode to be\ngenerally present in such heterojunctions, and that faster delocalization would lead to a reduced\nlifetime (i.e., less visible oscillation cycles). However, no oscillaitons are observed at all, also on\nshort time scales. Could it simply be that the mode is below detection capabilities?\nAuthors' answer. As Roelcke et al. argue, \"It is plausible that the atomically strong,\npredominantly out-of-plane THz electric near field efficiently couples to this drum mode\n(DM) by Coulomb interaction with charges in the heterostructure (for example, charge\ntransfer between Au and WSe ) […].\". This effect relies on the difference in workfunction\n2\nbetween the WSe and substrate, which will also affect the frequency of the drum mode.\norder of 1 eV, which can reduce the coupling to and frequency of the drum mode. Due\nto a comparable detection sensitivity as compared to Roelcke et al. we conclude that the\ndifference in the substrate's workfunction may be large enough to sufficiently suppress\nany experimental signature of the drum mode in our experiments, especially at small\ntimescales (<5 ps) where THz-induced tailing oscillations are present.\nReviewer: 4\nReviewer's comment.\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nAuthors' answer. We appreciate you taking part in Nature Communication initiative to facilitate\ntraining in peer review."}]

Decision letter