Ultrafast plasmonic photoemission in the single-cycle and few-cycle regimes

Due to the highly increased interest in the development of state-of-the-art applications of photoemission in ultrafast electron microscopy, development of photocathodes and many more applications, a correct theoretical understanding of the underlying phenomena is needed. Within the framework of the single active electron approximation the most accurate results can be obtained by the direct solution of the time-dependent Schrödinger equation (TDSE). In this work, after a brief presentation of a numerically improved version of a mixed 1D-TDSE method, we investigated the characteristics of electron spectra obtained from the surface of metal nanoparticles irradiated with ultrashort laser pulses. During our investigation different decay lengths of the plasmonic-enhanced incident field in the vicinity of the metal were considered. Using the simulated spectra we managed to identify the behavior of the cutoff energy as a function of decay length in the strong-field, multiphoton and transition regimes.

www.nature.com/scientificreports/ Throughout this paper atomic units (at.u.) are used ( = m e = e = 1 ), and for the incident laser fields we used Gaussian pulses having the same central wavelength = 800 nm (IR pulse) but different peak intensities and pulse lengths.

Theoretical method and numerics
It was previously shown in preceding works 20 that photoelectrons induced by laser pulses from different plasmonic nanoparticles have the larger kinetic energies the closer to edges of the nanotargets they escaped from. This observation can be explained by the fact that the edges of a nanorods/nanotriangles, or the tips of a nanotip are acting like nanoantennas for the incident EM field, i.e., under certain conditions (due to plasmonic effects induced by the field itself in the metal) the incoming E(t) electric field in the close vicinity of the target can be enhanced. The resulting locally enhanced field can be given as E loc (z; t) = Q(z)E(t) with Q(z) > 1 being the enhancement factor, which decays exponentially as a function of the distance from the surface.
The dynamics of these high energy electrons that probe the enhanced electric fields in the vicinity of the nanotargets' sharp edges are modeled by assuming a simple 1D scheme. Within this model the motion of the photoelectron points in the 0z direction that is perpendicular to the local surface and it is induced by the electric component of the laser field parallel to 0z axis. Here, 0z represents the symmetry axis of the laser-nanotip/edge system that is assumed obeying cylindrical symmetry around the normal of the local surface.
The time-dependent Schrödinger equation (TDSE). By taking into account the aforementioned considerations and in order to study more appropriately the dynamics of the laser emitted electron from the investigated nanorod (or nanotriangle) we employed the one-dimensional time-dependent Schrödinger equation (1D-TDSE) for an active electron located on the surface of the metal: where T = − ∂ 2 z /2 is the kinetic energy operator of the electron represented by the �(z; t) wavefunction ( ∂ z,t = ∂/∂z, t).
We consider a gas of quasi free electrons, for which the boundary of the metal means a confining potential. Using single active electron approximation, the electron experiences the following potential ( Fig. 1a): where E F is the Fermi energy, W is the work function of the metal ( W Au ≃ 5.1 eV in our case for gold), while the parameter β represents the screening constant of the image potential, whose value was set to β = 0.6 as showed to reproduce numerical results in good agreement with the experimental data 21 . Considering a Gaussian incident laser pulse with the electric field component (Fig. 1b): the laser-electron interaction potential was given within the dipole approximation and using its length gauge form: V le (z; t) = � r � E(t) ≡ zE(t) . In Eq. (3) τ p = τ FWHM was the pulse duration (in the full-width at half of the maximum intensity), ω 0 the central circular frequency, and ϕ CEP represents the carrier-envelope phase of the pulse. Furthermore, we assumed that the shielding-effect of the surrounding surface electrons onto the studied electron is sufficiently large, hence, the effect of the incident E(t) field inside the metal can be neglected:  www.nature.com/scientificreports/ for z < 0 . The field-enhancing plasmonic effects were considered by multiplying the incoming field with the Q(z) decaying enhancement-curve, that fixes the laser-electron interaction term to V le (z; t) ∼ = zQ(z)E(t) = zE loc (z; t).
From the numerical point of view, we are aware of that a 3D-TDSE calculation would provide us a much more complete study of the investigated photoemission phenomena, for instance also including the possibility of obtaining the angular energy distribution of the ejected electrons; nonetheless, such investigations concerning laser pulses with long wavelengths ( ≥ 800 nm) and durations ( τ ≥ 2.6 fs) would imply unfeasible expensive computational demands. In contrast, by using the 1D procedure, the dominant photon-induced effects, including the photoemission, could be still properly investigated with a less time-consuming procedure.
The mixed split-operator and Crank-Nicolson (SOCN) approach. For the numerical representation of the wavefunction (WF) we employed the finite-difference grid method. Within this approach, the potential energy terms, V(z) and V le (z; t) , are represented by diagonal matrices, while the second order differentiation operator ∂ 2 z �(z) /�z 2 by a tridiagonal matrix. Thus, their sum, i.e., the Hamiltonian will be a sparse matrix having nonzero elements only on the three main diagonals. In the simplest approach the Û (t + δt, t) = exp{−iδtĤ(t)} time-evolution operator is temporally discretized ( δt → �t ) and approximated with its Taylor-expanded form Û ≃ exp{−i�tĤ(t)} = n max n=0 (−i�tĤ) n /n! . The main drawback of this approach consists in the fact that the stability of the WF's temporal propagation can be achieved only for large values of n max terms (a radiation field dependent parameter), which would imply a large number of matrix power calculations in each time-step, while the time propagation itself remains non-unitary by definition (i.e., the WF norm, , is not conserved). This non-unitarity can be avoided with the use of the Crank-Nicolson (CN) procedure 24 , according to which a forward [ �(t ′ ) = Û (t ′ , t)�(t) ] and a backward [ �(t ′ ) = Û † (t ′ , t + �t)�(t + �t) ] propagation step to an intermediate state at time t ′ = t + �t/2 is performed, and the WF is evaluated as In this way the unitarity is ensured. Although the CN is an unconditionally stable and unitary approach it still necessitates a considerable amount of computation time, because for each integration timestep the calculation of a matrix inverse is needed. This shortcoming becomes even more pronounced with the increase of the wavelength and/or intensity of the radiation field due to the increased active space the photon-emitted electron travels in.
In order to overcome these shortcomings (e.g., multiple matrix inversions) one may use some spectral method 22,23 , that on the other hand would imply additional convergence tests regarding the required number of basis functions for the correct dynamics; or, for instance, just a simple split-operator method, which however could not guarantee the stability during the temporal propagation. Hence, in order to avoid additional convergence tests, or multiple matrix inversions, and to ensure stability in the laser-induced dynamics, we considered the combination of the two methods: (i) the CN approach 24 and (ii) the split-operator 25 technique. According to this mixed approach (SOCN) 21 the WF is divided into a � 0 (z; t) field-free and a � (z; t) field-perturbed part: = 0 +˜ . The solution in time moment t i + t of the stationary case ( i∂ t � 0 =Ĥ 0 � 0 = [T +V (z)]� 0 ) within the CN approach read as: For the field-perturbed part the following first-oder non-homogeneous linear differential equation (DE) can be deduced which is generally solved by summing up a particular solution ( ˜ p ) and the general solution of the corresponding homogeneous DE, i.e., where it is assumed that during the short t time interval V le (t) ≃ V le (t + �t) ≈ V le (t + �t/2) . In our approach, that considers this assumption, we deduced the solution of Eq. (5), first, by applying the CN scheme of the full WF: Then, by equalling Eqs. (6) and (7), and using �t 2 = �t/2 along with the split technique the field-perturbed part in time step t i+1 = t i + t can be given as After employing the first-order Taylor expansions of the exponentials inside the square brackets, and of the exp{i t 2Ĥ0 } term (in lhs of Eq. (8)) one finally obtains In accordance with the CN scheme the field-free evolution operator is obtained by evaluating the Û CN 0 = (Û † 0 ) −1Û 0 matrix-matrix multiplication, while by adopting the split technique the external timedependent field ( V le ) is detached from the full Û CN operator. As a consequence, the inverse calculation of the evolution operator in each integration step is avoided, saving this way an appreciable amount of computation time. Using SOCN the Û CN 0 (�t j ) matrices should be calculated only once, and stored in the computer's memory at the very beginning of the simulation. Here �t j = �t init /2 j−1 ( j = 1, 2, . . . ) represents the used discrete time intervals in the adaptive temporal propagation scheme ( t 1 = 10 −3 at.u.).
The initial state's wavefunction and the solution of the TDSE. The initial state WF [ � 0 = �(z; t = 0) ] was calculated by the direct diagonalization of the field-free Ĥ 0 =T +V (z) Hamiltonian (large and sparse matrix) with the use of the Scalable Library for Eigenvalue Problem Computations (SLEPc) package 27,28 . Before initiating the TDSE runs (solving Eqs. (4), (10)), we determined the optimum value for the numerical discretization parameters: z (gridpoints separation) and ± z max (grid-size). This was achieved using the SLEPc, and the value of the two discretization parameters was considered optimum when the calculated initial state WF converged. The SLEPc package can deliver that WF whose eigenenergy is the nearest to a usergiven target eigenenergy. Since we considered for the active electron a surface electron, for the target eigenvalue the energy of E F + W was used. Once � 0 (z) corresponding to the input energy was obtained, its accuracy was tested by defining the convergence parameter where i ≥ 1 , and for the initial value of z 0 we arbitrarily chose 0.8 at.u., while for i > 0 values from �z i ∈ {0.7, 0.6, 0.5, . . . } . For the optimum gridpoint separation distance we fixed the value of z = 0.5 at.u. when the d � 0 conv (�z) < 10 −8 condition was fulfilled. In order to diminish the undesired spurious reflections from the 'walls' of the simulation box during the temporal propagation, and to monitor the extent part of the WF, we employed at edges of the simulation grid complex absorbing potentials (CAP)of the form 29 : where z cut is the first point of the absorbing region. The size of this region was z max − z cut = 50 at.u. A similar CAP was used in the z < 0 direction as well, however the absorption/reflection here is much less probable since the accelerating effect of the laser field inside the metal is not considered. By performing an initial testing for the case of the highest field strength (with peak intensity I 0 = 80 TW/cm 2 ) and longest pulse (3 optical cycles), we calculated the amount of the absorbed WF at these far distances, and noticed that by increasing the simulation space to z max = 1500 at.u. ( ∼ 79.3 nm) in both directions ( z < 0 and z > 0 ) the sum of the two cumulative absorbed WF norms stayed below the value of 10 −10 . This is an upper limit for all the considered laser intensities and pulse lengths, thus the full information of the system was maintained inside the simulated area. In this way a sufficiently large simulation 'box' of z ∈ [− 1500, 1500] at.u. was constructed ensuring that even for the longest and the highest intensity pulse the relevant part of the ionized WF did not reach the boundaries (at z = ± z max ). In this way nonphysical simulation events that could significantly alter the characteristics of the final WF was prevented. During the integration of the TDSE in each time step the value of t was initially set to 10 −3 and later on adaptively modified (decreased iteratively by half) as long as the temporal propagation error defined as exceeded the value of 10 −8 . However, since the calculation of (Û † 0 ) −1 is a strongly time-consuming part of the simulation, in order to save even more CPU time we tried to reduce the number of required numerical matrix inversions during the complete temporal propagation as much as possible. This was achieved by using for the time-stepping t = 1 2 × 10 −3 , a value for which the total number of required inversions was kept acceptably low even for the case of the highest peak intensity (i.e., when the induced electron dynamics was changing in the fastest way).

Photoelectron spectra and wavefunctions
By applying the aforementioned considerations in the solution of the TDSE, using the optimum values for the discretization parameters we solved Eqs. (4)-(10) for different peak intensity and central wavelength pulses in the single and few-cycle regime. As a result, the final electronic WFs � (z; t = τ ) and their associated photoelectron spectra (PhES) were calculated. The PhES were obtained by projecting the ionized (departed) part of the final www.nature.com/scientificreports/ wavefunction � + (τ ) =�(z > 0; τ ) onto the continuum states |ψ k (z)� = (2π) −1/2 e ikz (free plane waves) with kinetic momentum k: where W k (τ ) represents the ionization probability density at the end of the laser field.
In the first step, we omitted the decaying amplitude of the plasmonic field enhancements and by fixing the field strengths through the whole coordinate space to E(z; t) = E 0 (t) we oriented our investigation onto the final dynamics obtained in the proximity of three different intensity regions: multiphoton γ > 1 ( I 0 = 5 TW/cm 2 ; γ = 2.8 ), transition 30 γ ∼ 1 ( I 0 = 40 TW/cm 2 ; γ = 1.02 ), and strong-field γ < 1 ( I 0 = 80 TW/cm 2 ; γ = 0.72 ) regime. We plotted in Fig. 2 the square of the final WF's absolute value (thicker red curves) and the PhES (thinner black curves) calculated according to Eq. (14) in this three distinguished regimes for different pulse lengths: n OC = 1 , n OC = 2 , n OC = 3 optical cycles at the FWHM of the intensity ( n OC = 1 corresponds to 2.66 fs).
As one may observe by increasing the incident field's peak intensity, a plateau-like region of electrons start to build up in the PhES in the case of each pulse lengths. As long as for γ = 2.8 the shape of the obtained spectra are more triangle-like, when the value of the Keldysh-parameter drops (intensity is increased) higher and higher energy (momenta) electrons start to appear with more or less the same probability in the middle region of the momentum maps: 1 < k < 1.5 at.u. for γ = 1.02 Fig. 2d-f, and 1.2 < k < 2 at.u. for γ = 0.72 Fig. 2g-i. The data show that these values do not depend on the pulse duration but exclusively on the strength of the driving field (observe the similar plateau for each intensity and different n OC ). Moreover, as expected, the start and endpoint of this regions correspond to 2U p and 10U p energy values, respectively. U p = E 2 loc,0 /4ω 2 0 is the ponderomotive (cycle-averaged quiver) energy of the photoelectron in the laser field, where ω 0 = 2πc/ represents the circular www.nature.com/scientificreports/ frequency and c is the speed of light. This means that upon 2U p the directly ejected and accelerated electrons are present (predominantly for the case of γ = 2.8 ), but by further increasing the pulse intensity higher energy electrons emerge in the PhES whose energy is greater than 2U p . This translates to the fact that a considerable portion of the ejected electrons are re-driven and accelerated back by the oscillating field onto the metal's surface, where they collide with this by experiencing a field-modified potential. As the external field becomes reversed, the returning electron wave packet-provided that the potential-barrier's distortion is strong enough-may experience a broaden and higher potential barrier (compared to the one when it was released in the continuum), and will undergo a forced and rapid deceleration phase. The newly emerged outward force due to the built-up barrier adds to the force exerted by the momentary field of the recently arrived optical half-cycle, and according to this the photoelectron will be accelerated even more in the outward direction due to the newly developed and increasing net force. In this way the electron may gain as much as ∼ 10U p kinetic energy at the end of the optical half-cycle (and eventually at the end of the laser pulse itself). The presented data in Fig. 2 undoubtedly proves the presence of this secondary dynamics of electrons starting from γ ≈ 1 and below 1. As seen, the increase of the pulse length does not have an important effect on this phenomenon, however since the more optical cycles are involved the more the electrons are ejected, and due to the increased number of electron wave packet trajectories interferences of these may recurrently occur bringing additional small wave-like features on top of the PhES. While the cutoff energy is the same regardless of the pulse length (as already discussed previously), with the increase in the number of oscillating optical cycles the spectra show more and more sub-features, subcomponents. These presumably correspond to above-threshold-photoemission (ATP) peaks, and one can see that for higher incident intensities, when the plateau starts to build-up in the PhES, this emerging plateau region is basically represented by a train of ATP peaks having about the same probability amplitude. The appearance of these wiggle-like features is even more straightforward by looking at the curves of the WFs: with the increase in the number of applied optical cycles the details are more obviously present. In addition to that, if one considers the shape of the WFs may note that for n OC = 1 well distinguishable features/peaks can be identified, whereas by increasing n OC next to the increased extent of the electron wave packets its higher complexity can also be observed.
Another observation is, that until a pulse is shorter and has a smaller number of ionizing optical cycles, the shape of the WF and PhES are more similar to each other, both having well distinguishable segments in their curves. On the other hand, the longer pulses will bring more discrepancies into this correspondence, mostly because the characteristics of the WF's shape is slightly more affected by this than the shape of its momentum counterpart (i.e., the PhES), where by the integration over the space may smooth out some of the wiggle-like features present in the curves of the electrons wave packets. Nevertheless, since a fine correspondence between the shape of the ionized WF and PhES can be determined for ultrashort (single, few cycle) pulses, this could promise novel possibilities for evaluating and comparing experimentally acquired PhES by simply calculating and investigating the electron's WF.

Plasmonic effects and the cutoff law
In quite a many atomic physics paper a well known formula is used for defining the cutoff energy of the laserionized photoelectron's, which reads as: where next to the 10U p (with U p ∼ I 0 2 ) a 0.537 W quantum correction term appears with W being the ionization energy of the atom. Contrary to photoemission from metals, this formula was deduced for atomic potential and external radiation field that was considered homogeneous throughout the space.
Certainly, in the present work the potential (Fig. 1a) is different from the atomic one and in the case when plasmonic effects in the vicinity of the target are also taken into account the radiation field loses its homogeneity as well. In order to check the applicability of the referenced formula (Eq. (15)) for the plasmonic metallic targets we performed several runs where the decaying profile of the plasmon-enhanced field was considered in the three ionization regimes (multiphoton, transition, and strong-field). Here we considered an exponential decay along the 0z direction of the local E 0,loc field by multiplying it with the exp{−z/l f } factor, where l f represents the decay-length of the radiation field. The δ = l f /l q was introduced to quantitatively indicate the relation of the decay-length's value to the amplitude of the quiver motion maintained by the laser field, l q = qE 0 /mω 2 .
In Fig. 3a-c we plotted the final spectra obtained for a set of the field-decay parameter values: δ ∈ {0.2, 0.5, 1, 5, 10, 50, 100} ; and also for the case when the laser-field of n OC = 1 was considered homogeneous along the 0z direction at a constant value of E(z; t) = E(z = 0; t) = E(t) , which corresponds to δ = ∞ . As one can see, irrespective of the local value of γ in the PhES obtained for δ ≤ 1 the ionization yields are low and can be characterized more with a triangle-like shape, i.e., given by the presence of the so-called roll-off electrons. Even in the case of the highest intensity where a large-scale plateau is already expected, only a short region of lower energy electrons are showing up with similar probability amplitudes. It is also obvious that as the value of δ becomes higher (the decay-length longer) the ionization probability density curves begin to rise until a final level is approached starting from δ ≥ 5 . This observation can be noticed more clearly from Fig. 3d where the calculated ionization yields (the integral of dP/dk over k is performed) converge to the value obtained for δ = ∞ homogeneous case. This convergence starts to take effect above δ = 1 , and is slightly faster for smaller γ values (higher intensities). By calculating the cutoff points given by Eq. (15) in each separate regime, we see that the expected cutoff momenta are attained only for δ > 10 , which proves that the cutoff formula derived for atomic targets and spatially homogeneous oscillating fields is valid only above δ ≥ 10 , and can be safely used starting from δ ≥ 50 for 800 nm ultrashort IR pulses. www.nature.com/scientificreports/ Since for the lower δ values the correct identification of the cutoff point is hard or even impossible to achieve, we compared the obtained spectra when δ = ∞ to the case of the homogeneous field ( δ = ∞ ) by defining the following equation: which quantity tends to 1 for δ → 0 and to 0 for δ → ∞ . As one can see in Fig. 3e for δ ≥ 5 this quantity for each intensity cases is already below 0.1, however in subfigures (a-c) we can observe that the cutoff point is still relatively far from the 10U p value. Only starting from δ ≥ 50 we see very similar final spectra and cutoff point for each regimes. Moreover, we can observe that for higher local field amplitudes ( γ loc = 0.72 ) the cutoff for δ = 50 is located further away from the 10U p value than in the case of lower intensity ( γ loc = 2.86 ). This observation indicates that the decaying profile of the plasmon enhanced electric field has more impact on the final spectra when the incident field has higher amplitude.
The explanation of this finding is encoded in the relationship between the extent of the emitted electronic wave packet and the decay-length of the plasmonic field. For lower δ values the incident field's effect gets suppressed for larger distances measured from the nanotarget surface, and the wave packets ionized in the first part of the pulses could reach far regions from where cannot be driven back by a fast decaying driving field, ruling out this way the recollision. Moreover, a more abruptly decaying field would transfer less energy for the motion of the electron, than a homogeneous constant field would do through the whole coordinate space. Another aspect constitutes in the lower yields obtained for the same local intensity but different decay-lengths. The reason why one gets smaller amount of ejected electrons after the completion of the laser pulse can be easily explained considering Fig. 3f, which shows the distortion suffered by the potential barrier at the top of the main/central optical half cycle. It is obvious that for lower δ parameter values the width of the field-distorted barrier is larger, hence the probability of ionizing an electron from the initial state drops. This translates to fewer electron-birth events in the continuum, decreasing in accordance with this the final ionization yields. However, in real-life scenarios  26 . In (d) the total ionization rates obtained for different Keldysh and delta paramater values. In (e) the relative differences between the spectra calculated for a certain δ < ∞ to the PhES obtained when the field decay was absent ( δ = ∞ ) for the case of the three ionization regimes. Subfigure (f) illustrates the shape of the electron's V(z) potential distorted by the maximum of the laser field for the cases of different δ decay-parameter values. www.nature.com/scientificreports/ it may happen that due to the field-disturbed surface even for the case of lower δ values the active electrons could reach levels of higher excited states before the ionization takes place. From this higher energy levels the probability amplitude of escaping a broader barrier is increased and the same the electron yields in the continuum. This may alter a bit the real convergence seen in Fig. 3d, however the position of the cutoff region wouldn't be impacted much, since it is more affected by the driving field's strength in the continuum, thus our statement regarding its dependence upon the value of δ can be considered to be still valid.

Conclusions and outlook
In the first part of this work, we presented a computationally efficient method to study the laser-initiated dynamics from metallic nanotarget surfaces for a single active electron. The method was based on the direct solution of the one-dimensional Schrödinger equation represented on a finite-difference grid. We showed that by the diagonalization of the field-free Hamiltonian-matrix a set of (initial) eigenstates can be obtained at once at the start of our investigation. By considering the active electron initially on the Fermi-level we calculated the final wavefunction and ionization spectra resulting from the interaction with external (800 nm) laser fields of different durations and peak intensities, by covering the three intensity regimes: multiphoton, strong-field, and transition. It was shown that by increasing the incident field intensities and by entering more in the strong-field regime (with the Keldysh-parameter γ ≤ 1 ), a plateau of higher energy electrons started to appear in the calculated spectra (characterized by a traceable cutoff point, starting from which the highly energetic electrons start to disappear from the spectra). These are strong-field effects that results in a detectable amount of high energy electrons in the measurements that also encode information regarding the enhanced local fields in the close vicinity of plasmonic nanotargets. In the final part of the paper, we showed the applicability of a cutoff-law that was deduced for atomic potentials also for the case of photoemission taking place from plasmonic nanoobjects considering different field enhancement profiles (decay lengths). We showed that the referenced cutoff-law can be safely utilized for photoemission when the field-enhanced local fields have not so abrupt decrease ( δ ≥ 50 ) in the near vicinity of the surface.