a, Terahertz voltage waveforms used in the simulations. We approximate the time trace of the THz voltage with the THz electric-field waveform measured in the far field by electro-optic sampling. In all plots, the pink curve corresponds to a waveform with a peak of −2.05 V and the navy-blue curve corresponds to a waveform with a peak of −2.22 V. The scaling factor and resulting peak voltages are determined by fitting the shape of the experimental autocorrelation (Fig. 3a) as follows. We simulate the autocorrelation by taking the sum of two THz pulses and calculating the total rectified current from the resulting transient as a function of the delay time between the pulses. b, The dI/dV characteristic needed for the simulation is obtained by modelling the HOMO and LUMO molecular-resonance peaks in the experimental dI/dV curve (Fig. 2b, centre) with two Gaussians. This dI/dV curve is then scaled up to account for the significantly lower tip height for THz-STM imaging. Note that the scaling factor does not affect the shape of the autocorrelation. Blue line, centre of HOMO Gaussian, −1.93 V (0.62 V full width at half-maximum, FWHM). Red line, centre of LUMO Gaussian, 2.06 V (0.58 V FWHM). c, Resulting current response induced by the THz voltage waveform in a when applied to a junction with a dI/dV relation defined by b. The asymmetry of the THz voltage pulse leads to a much larger current response in the negative bias direction than in the positive bias direction. Inset, induced current response during the negative crest of the THz voltage waveform. Pink curve, 120-fs FWHM; navy-blue curve, 140-fs FWHM. d, Total number of electrons that have tunnelled across the junction, calculated from the integral of c. Negative numbers refer to the negative bias direction, that is, electron tunnelling from the HOMO to the tip. Inset, rise of the rectified electron signal during the negative crest of the THz voltage waveform. The rise time from 10% to 90% of the maximum signal is 115 fs for the pink curve and 130 fs for the navy-blue curve. The simulated autocorrelation for the −2.22 V peak (navy-blue curve) provides the best fit to the shape of the measured autocorrelation. The dI/dV curve is scaled such that the autocorrelation peak at −2.22 V matches the measured peak of approximately −0.75 electrons per THz pulse. The peak voltage of the THz pulses used in the imaging configuration is then determined by finding the THz peak for which the simulations yield this number of rectified electrons per THz pulse. For example, in Fig. 2e the maximum observed signal is approximately −0.58 electrons per THz pulse, and the best agreement is found for −2.05 V (pink curve). We note that the simulations are based on the assumption that the THz pulse modulates the bias of the junction quasi-instantaneously. In our simulations we disregard any blocking of tunnelling or other effects resulting from the finite time until an electron from the substrate refills the molecular state. This blocking is expected to become important if the THz current reaches or even exceeds one electron per pulse.