Third step: calculation of electron trajectories

In the third step, the electron trajectory simulation was performed. The total yield is obtained by integrating the calculated emission current over energy, time and space. The number of electrons per pulse was 850 at the highest laser intensity. The initial electron emission times, energies, positions and directions were determined based on the Monte Carlo method by using the calculated emission current distribution. The tunnelling emission was delayed by the Keldysh time, τ_K, multiplied by a factor C_τ. Here, τ_K is expressed as τ_K = γT₀/4π where γ is the Keldysh parameter and T₀ is the laser oscillation period⁵⁷. The Keldysh time varies with position and time of electron emission because it depends on the instantaneous laser fields and the effective barrier height reduced by DC fields. With this model, we effectively determined the start position and time after tunnelling dynamics. Whether the emitted electrons can be driven back to the surface by the oscillating laser field depends on the tunnelling delay time in our model. C_τ was adjusted in such a way that the largest number of electrons is redirected to the surface with the highest energy. This condition resulted in a value of 1.4 for C_τ for the case of Model C. If C_τ is set to zero then few electrons are re-directed towards the surface. In addition, the start velocities of the tunnelling electrons are calculated from the electron velocities inside the metal; the perpendicular component of the velocity with respect to the surface is set to zero when the electrons are at the tunnelling exit, and additionally we added a velocity change due to the vector potential experienced during the delay time. It should be noted that the time-dependent Schrödinger equation should be employed for accurate theoretical description of optical tunneling time as discussed in refs ⁹ and ¹⁰.

After determining the initial conditions, the electrons were propagated through the vacuum, where they experience four kinds of forces: 1. laser fields, 2. DC fields, 3. Coulomb forces between electrons (space charge effects), and 4. image charge forces from the tip surface. For calculating DC fields, we used the OpenMaXwell static solver with fictitious charges along the tip axis²⁸, which is basically the same technique as explained in ref. ⁵⁸. The tip is biased at approximately −3300 V and there is a grounded counter electrode 13 mm away from the tip apex, which represents the conditions of the experimental setup. The DC fields are set to 2.5 V/nm at the very top of the tip apex, which is a reasonable value when tunnelling emission due to the DC fields is observed. For simulating dynamics of image charges in the tip, we simulated the positions and amounts of classical image charges for a metallic sphere. The simplified spherical approximation of the tip is acceptable because image forces are important only very close to the tip apex. Tracking electron trajectories was done with including those image charges in the computational fields, and calculated Coulomb forces between all the electrons and image charges.

During propagation, the time steps were dynamically determined in such a way that the mean energy change of all the electrons for each time step was a constant. This constant determines the accuracy of the simulation; the cumulative absolute energy error is estimated to be about ±0.2 eV. The tracking-electron simulation stops when all the emitted electrons reach the counter electrode (pinhole plate). This entire process was repeated until enough statistics was obtained.

The electrons which are re-scattered from the tip experience the delay processes as illustrated in Fig. 1(b). In electron travelling in the metal, we assume a simple model where the electrons travel a distance of the order of twice the inelastic mean free path before reappearing in the vacuum. Here we applied an effective inelastic scattering rate when electrons have re-entered the surface. This rate is a global fitting parameter in our simulation, i.e. one common value for all laser intensities, which is adjusted in order to best describe the peak-to-plateau intensity ratios from the experimental data. The physics behind this quantity is complicated and includes contributions from inelastic forward and backscattering in the solid and also creation of secondary electrons^31,59,60. We obtained a value of 0.45 for best fit, and the number of electrons reappearing in the vacuum is effectively determined simply by using this inelastic scattering rate. Note that the elastic scattering probability is around 0.1–0.15 surface in our electron energy regime⁶⁰. Due to the concomitant creation of secondary electrons via cascading processes, the effective inelastic scattering rate at low energies can be higher than the elastic scattering rate, and a value of 0.45 appears still physically sound. Regarding elastic scattering, the typical time scale for travelling in the metal is very short, of the order of 1 fs or less for tungsten in our energy regime²⁷. Hence elastically re-scattered electrons can be considered as prompt, they contribute to the plateau feature or its high-energy tail. As discussed in the text, we neglected these processes.

The delay times due to the inelastic re-scattering processes were calculated by dividing the electron travelling length by the group velocity, taken from a free-electron model and using the effective mass of tungsten, for the round trip from and back to the surface. To determine the path length which electrons travel before the inelastic scattering takes place, we assume that the scattering events are stochastic. The probability distribution follows the Poisson distribution with its maximum positioned at the inelastic mean free path (IMFP) of tungsten⁶¹. The IMFP is obtained from ref. ²⁷. The thus calculated delay time is multiplied with a factor C_l to have an evolution of the peak feature similar to the observations. The resulting C_l was 1.8. This high value is considered to be due to the particularly low group velocity near energies of 5 eV above the Fermi level as calculated from dispersion relations²⁷.