Ultrafast strong-field photoelectron emission from biased metal surfaces: exact solution to time-dependent Schrödinger Equation

Laser-driven ultrafast electron emission offers the possibility of manipulation and control of coherent electron motion in ultrashort spatiotemporal scales. Here, an analytical solution is constructed for the highly nonlinear electron emission from a dc biased metal surface illuminated by a single frequency laser, by solving the time-dependent Schrödinger equation exactly. The solution is valid for arbitrary combinations of dc electric field, laser electric field, laser frequency, metal work function and Fermi level. Various emission mechanisms, such as multiphoton absorption or emission, optical or dc field emission, are all included in this single formulation. The transition between different emission processes is analyzed in detail. The time-dependent emission current reveals that intense current modulation may be possible even with a low intensity laser, by merely increasing the applied dc bias. The results provide insights into the electron pulse generation and manipulation for many novel applications based on ultrafast laser-induced electron emission.

, E F and W are the Fermi energy and the work function of the metal respectively, and e is the electron charge. In eq. (1), the metal is approximated by jellium model, while the external fields are assumed cut off abruptly at the surface. The "sudden screening" of external fields may be justified, because the field penetration depth (i.e. skin depth) δ m = (2/ωμ 0 σ) 1/2 of metal is typically much smaller than the laser wavelength (e. g. for gold, the conductivity σ = 4.1 × 10 7 S/m, for 800 nm laser, δ m = 4.06 nm). To calculate the emission probability, we solve the time-dependent Schrödinger equation, where ψ(x,t) is the complex electron wave function, m is the electron mass, and the potential Φ (x,t) is given in eq. (1). Following the coordinates transformation given by Truscott 33,34 , it is found that eq. (2) has an exact solution for x ≥ 0 (see Methods), where ε is the initial energy of an incident electron, T n is a constant representing the transmission amplitude, and G n (x, t) is given by, 2 . It can be easily verified that eq. (3) is a solution of eq. (2) by direct substitution. Equation (4) represents an outgoing wave traveling to the right, with Ai and Bi being the Airy function of the first kind and second kind, respectively. It is clear that both G n (t) and the second exponential term in eq. (3) are time-periodic functions with period of 2π/ω, thus, eq. (3) represents the superposition of electron plane waves with energies of ladder states with spacing of the photon energy ω  coupled to the initial energy at ε, as shown in Fig. 1b. The transfer of incident electrons at ε to sub-bands of  ε ω + n is made possible by processes such as multiphoton The applied dc electric field is F 0 , the laser electric field is F 1 cosωt. The initial electron energy is at ε. (b) The resulted ladder states for electrons with spacing of the photon energy ω  coupled to the initial electron energy at ε. The incident electron wave is ψ i , the transmitted and reflected electron waves are ψ t and ψ r respectively, both are superposition of plane waves with energies at the ladder states.
Scientific RepoRts | 6:19894 | DOI: 10.1038/srep19894 absorption (n > 0) or multiphoton emission (n < 0). Note that these ladder eigenstates are similar to Volkov states for a time-periodic potential 4,20,35 , whose eigenenergy ε ω + n consists of a drift kinetic energy E n and a ponderomotive energy U p due to electron quiver motion, both of which are defined in the sentence following eq. (4).
For x < 0, the solution to eq. (2) may be written as, which is the superposition of an incident plane wave and a set of reflected waves with reflection amplitude R n , The reflected waves consist of an elastic reflected wave (n = 0) and waves with ladder eigenenergies representing multiple photon energy (n ≠ 0), as shown in Fig. 1b. Although the general solution in eq. (5) includes the reflected waves with n < 0, whose energies are below the initial electron energy ε (typically assumed to be the Fermi energy), it is verified that almost all the reflected current is through the initial energy level, i.e. n = 0. When  ε ω ≤ − / n , k n becomes imaginary and the corresponding reflected wave becomes evanescent, resulting in zero reflected current through the channels with ε ω ≤ − / n  . Note that in eq. (5) a plane wave is assumed for the incident electron, which turns out to be a good approximation to the bound-state electron wavefunction inside the metal solid, since its dimension in the longitudinal direction (perpendicular to the surface) is typically much larger than the electron wavelength. In fact, the laser field interacts only with the evanescent part of the electron states near the metal surface, which is found to be very similar to the bound-state wave function in the atomic case 29 .
The solutions in eqs (3) and (5) are matched at the metal-vacuum interface from the conditions that both ψ and dψ/dx are continuous at x = 0. By taking Fourier transform, we obtain, in nondimensional quantities, ε ε n n n n l n n l where δ is the Dirac delta function and, n n n , and a prime denotes derivative with respect to the argument. In eq. (7b), p n and z n represent the phase factor of the wavefunction in the nth state and of its spatial derivative at = x 0, respectively. In eq. (7a), P nl and Z nl are the lth Fourier coefficients of p n and z n in eq. (7b), respectively. The summation rule in eq. (6) physically means the conservation of probability at the metal-vacuum interface, = x 0. It is derived from the conditions that both the wave function and its first derivative be continuous there. The transmission amplitudes T n is calculated from eq. (6). The normalized emission current density is defined as w(ε, x, t) ≡ J t (ε, x, t)/J i (ε), which is the ratio of the transmitted current density over the incident current density, calculated by using the probability current density l . The emission current density is found in the normalized form as, n l i l n t n l nl nl nl l n n l n l n l 1 0 where w n describes the emission current density through the nth channel, with emitted electrons of energy ε ω + n , due to the n-photon contribution. Once the dc electric field F 0 , laser field F 1 , laser frequency ω, and metal properties (E F and W) are specified, the time-dependent and time-averaged emission current density may be calculated from eqs (8) and (9) respectively, with T n given by eq. (6), for any given initial electron energy ε. For simplicity, we take the energy of the initial electrons to be ε = E F , which is justified by the well-known fact that the probability of electron tunneling drops rapidly when electron energy decreases, with most of the electrons being emitted from sources near the Fermi level 20,23,31,32,36 . Numerical calculations have verified that the emission current drops exponentially as the initial electron energy decreases below the Fermi level.
Main results and discussion. Figure 2 shows the time-averaged normalized emission current density w n through the nth channel, under various combination of dc electric field F 0 and laser field F 1 , calculated from eq. (9). The laser wavelength was chosen to be 800 nm, or a photon energy of ω = 1.55 eV. The metal is assumed to be gold, with Fermi energy E F = 5.53 eV and the work function W = 5.1 eV. Unless stated otherwise, these would be the default values for the calculations in this paper. It is clear that the total emission current density w increases when either F 0 or F 1 increases. When both F 0 and F 1 are small (Fig. 2a), the dominant emission process is the three-photon absorption, which is consistent with the ratio of the work function to the photon energy ω / W  = 3.29, where the electrons at the Fermi level need to absorb at least three photons to substantially reduce the effective potential barrier W (Fig. 1a). As F 1 increases, higher order channels become important (Fig. 2b,c). In Fig. 2c, all the emission processes with 3 ≤ n < 10 contribute significantly to the emitted current. The energy distribution of emission current (i.e. the envelope in Fig. 2a-c), exhibiting a transition from narrow-peaked spectrum to a broad plateau-like spectrum, shows striking resemblance to the experimentally measured kinetic energy spectra 17 , for the transition from multiphoton regime to the strong field regime (c.f. the measured spectra in the insets of Fig. 2a in ref. 17). In Fig. 2c, the emission process with highest probability shifts to four-photon absorption, this is consistent with the well-known feature of channel closing when considering the laser field only (i.e. F 0 = 0) 28,29 . The minimum value of n increases if the given laser field F 1 increases, as determined by the drift kinetic energy of the electrons p , which has to be nonnegative. However, complete channel closing is not observed, especially for cases with larger dc electric field F 0 (Fig. 2d-i), where substantially more channels are opened for electron emission. The number of possible emission channels increases with both F 0 and F 1 . Besides multiphoton absorption (n > 0), electrons may also emit through processes such as optical direct tunneling (n = 0) and multiphoton emission (n < 0). As the dc field increases, the effective potential barrier becomes narrower (Fig. 1a), which would not only increase the probability for electron emission through existing channels, but also facilitate electron emission through more channels. For a given F 0 , increasing F 1 tends to shift the dominant emission process to channels with larger n, which is consistent with the trend shown in the experiment by Schenk et al. 37 . In particular, Schenk et al. 37 observed the electron emission yield at the S = 1 peak exceeds that of S = 0 peak at a laser intensity of ~2 × 10 11 W/cm 2 for the illumination of a sharp tungsten tip, where is the above-threshold order. Similar peak shift is noted in Fig. 2b-c, and in Fig. 2d-e. For a given F 1 , increasing F 0 tends to bring the dominant emission process closer to optical tunneling n = 0. Figure 3a shows the time-averaged normalized total emission current density w as a function of laser electric field F 1 (and laser intensity I) for various dc electric fields F 0 . When F 1 < 2 V/nm (I < 5.32 × 10 11 W/cm 2 ), the emission current density is well scaled as ∝ w F n n 1 2 , indicating the dominant n-photon process. The value of n decreases as the dc electric fields F 0 increases, which is consistent with the results in Fig. 2.
For the special case of zero dc electric field F 0 = 0, the solution in the previous section is simplified by replacing eq. (4) with  ξ , so that T n is still calculated from eq. (6), with eq. (7a) unchanged, but eq. (7b) replaced by . The normalized emission current density in eqs (8) and (9) is modified to read, Im sin 10 Equation (11) is identical to eq. (46) of Yalunin et al. 20 . The emission current density for F 0 = 0 calculated from eq. (11) is also shown in Fig. 3a, which serves as the lower limit for F 0 → 0. When F 0 = 0, four-photon absorption is the with n = 4. The abrupt drop of w around F 1 = 11.6 V/nm is due to the channel closing effect 20 , which is accurately predicted by letting , giving F 1 = 11.74 V/nm.
When F 0 is non zero but small (1 V/nm), the sharp peak of w near the channel closing region is smeared, however, there is still a clear transition from multiphoton absorption to optical tunneling emission, indicated by a slope change around F 1 ~ 10 V/nm (I ~ 1.33 × 10 13 W/cm 2 ). The corresponding Keldysh parameter, which is usually used to characterize this transition, is γ ω = / = / W U F 2 p 1 ≈ 1.8. For F 1 < 10 V/nm, the emission current is well scaled as ∝ w F n 1 2 with n = 3, indicating the dominant three-photon absorption. For F 1 > 10 V/nm, the scale ∝ w F n 1 2 no longer applies, as the emission enters the optical field emission regime. When F 0 becomes large (≥ 2 V/nm), the slope of w changes gradually as F 1 increases, indicating that multiple processes of high photon numbers contribute, and that the dominant process changes as F 1 increases, consistent with Fig. 2. The insets in Fig. 3a shows the contribution of different processes for the case F 0 = 5 V/nm. The dominant process changes from single-photon absorption (n = 1) when F 1 < 2 V/nm (I < 5.32 × 10 11 W/cm 2 ) to two-photon absorption (n = 2) when F 1 > 2 V/nm (I > 5.32 × 10 11 W/cm 2 ). The individual n-photon absorption follows closely the scaling ∝ w F n 1 2 for n ≥ 1, when F 1 is small. The small dips along F 1 are due to the channel closing effects 20,28 . The direct tunneling process (n = 0) is almost independent of F 1 . The single-photon emission process (n = − 1) follows the same scaling as the single-photon absorption (n = 1), with ∝ w F n 1 2 , but with a substantially smaller emission current density. Figure 3b shows the total time-averaged emission current density w as a function of dc electric fields F 0 for various laser electric fields F 1 . When F 0 is small (< 2 V/nm), w is insensitive to F 0 for small F 1 (≤ 5 V/nm), because the dominant emission mechanism in this regime is multiphoton absorption, which presumably is independent of the dc field F 0 . For larger F 1 (≥ 7 V/nm), w increases with F 0 , since optical field emission would be important and applying a larger F 0 will further assist the tunneling by lowering the potential barrier. When F 0 is large (> 3 V/nm), Fowler-Nordheim like field emission 31 from the dc electric field becomes important, thus w increases with F 0 , for all the values of F 1 . As F 0 further increases (> 8 V/nm), w approaches that of field emission in static electric field (i.e. F 1 = 0), which, when calculated from eq. (9), gives identical emission current density to that of the known exact solution 31,32 . This may serve as one validation of our analytical solution. Note that eq. (12) reduces to the Fowler-Nordheim (FN) equation, , which is also plotted in Fig. 3b. In the regime of large F 0 , a larger F 1 increases <w> even further. Note that the trends of w in Fig. 3b are very similar to those in the voltage-and laser power-dependent electron flux measured experimentally by Ropers et al. (c.f. Fig. 2c of ref. 2). The insets in Fig. 3b shows the contribution of different emission processes for the case F 1 = 3 V/nm. The dips for the individual n-photon absorption is due to the channel closing effects 20,28 . As F 0 increases, the dominant emission process changes continuously, from three-photon absorption (F 0 < 2 V/nm), to two-photon absorption (2 V/nm < F 0 < 5 V/nm), to single-photon absorption (5V/ nm < F 0 < 10 V/nm), and eventually to direct tunneling with n = 0 (F 0 > 10 V/nm). Figures 4 and 5 show the time-dependent total electron emission current density w(t), under various combination of dc electric field F 0 and laser field F 1 , calculated from eq. (8). The average emission current is the same as those shown in Fig. 2. The emission current is periodic in time, producing electron bunches with the same period as the laser field 2π/ω. This is consistent with the experimental observation of periodic variation (at the laser period) of the electron yield due to the carrier-envelope phase (CEP) effects 17 . It is interesting to note from eq. (8) that the emission current from each individual channel, However, this is not the case for the total emission current, as shown in Figs 4 and 5. The complicated shape and phase shift of the total emission current are due to the highly nonlinear cross-terms between different channels (i.e. terms with n ≠ l in eq. (8)). Note that the instantaneous emission probability does not necessarily go to zero when the net electric field is pointing away from the surface, because of highly nonlinear dynamics of electrons near the metal-vacuum interface. This nonlinear process has been thoroughly discussed in previous work (e.g. Figs 5 and 6 of ref. 20, or Fig. 5 of ref. 38). In general, increasing the laser field F 1 would increase the peak-to-average ratio of the emitted current, meanwhile reducing the FWHM of the current pulse. In contrast, increasing F 0 tends to increase the FWHM of the current pulse. Most importantly, for a given laser field F 1 , increasing the dc electric field could substantially lower the valley of the current pulse, sometimes even close to zero (Figs 4d,h,i and 5d,h,i), thus giving an intensive modulation of the emitted current. When the dc bias F 0 is sufficiently high, profound modulation of the current emission is even possible for a relative small F 1 , (Figs 4d,g  and 5d,g), i.e. for a low intensity laser. It is interesting to note that, at locations x larger than 10 (beyond the localized surface current oscillations, as shown in Fig. 4), the current keeps approximately a fixed temporal profile with only a phase shift, especially when dc bias F 0 is larger than laser field F 1 (Figs 5a,d,g,h). Physically, the phase shift is primarily due to the drift and acceleration motion of the electrons under the influence of dc and laser fields. The time dependent emission current and its relative phase shift would be important for the study of carrier envelope phases (CEPs) 1,17,19,39 . Figure 6 shows the effects of laser frequency (or wavelength) on the emission current density. When the dc electric field F 0 is small, the effects of laser frequency become more pronounced. Figure 6a plots the current emission w as a function of laser field F 1 (and laser intensity I), for various combinations of laser wavelength and F 0 . When F 0 = 1 V/nm, the emitted current follows closely the power scale ∝ w F n 1 2 in the multiphoton regime when F 1 is small, where n is determined by the integer round off of the ratio ω / W  = 2.06, 3.29, and 4.11, for wavelength of 500 nm, 800 nm and 1000 nm, respectively. The laser field F 1 for the transition from multiphoton absorption to optical field emission increases with the laser frequency, consistent with the Keldysh parameter, γ ω = /F 1 . For larger F 0 , the effects of laser frequency become less important. When F 0 = 10 V/nm, the electron emission depends little on the laser frequency. Figure 6b confirms that the effects of laser frequency are important only when the dc bias is small. Figure 6c plots the current emission w as a function of laser frequency, for various combinations of F 0 and F 1 . When F 0 = 1 V/nm, there are distinct peaks near ω = 1/n (n > 1), corresponding to the n-photon process, with decreasing electron emission current as n increases. For ω > 1, the energy of single photon exceeds the potential barrier W, thus the electron emission process is similar to that of over-barrier ionization, where the emission current density decreases as frequency increases. When F 0 = 5 V/nm, the n-photon peaks are smeared for ω < 1; however, the emission current densities for ω > 1approach those of F 0 = 1 V/nm, indicating the over-barrier absorption. When F 0 = 10 V/nm, the emission current density becomes almost independent of the laser frequency for ω < 1, since the dominant emission process is the Fowler-Nordheim like field emission. When ω > 1, both dc field emission and over-barrier emission would be important, where the dominant process is determined by the strength of the laser field F 1 . It is important to note that the maximum emission current density is always around ω = 1. Figure 7 shows the effects of the barrier height (or work function) W on the emission current density, for E F = 5.53 eV and ħω = 1.55 eV (λ = 800 nm). In general, the emission current density increases as the work function decreases. When both F 0 and F 1 are small, there are peaks around W/ħω = n, corresponding to the resonant n-photon absorption. Similar results are obtained by previous numerical calculation 24 . When F 1 increases, the oscillating features persist, however, the location of the peaks shifts, due to both channel closing effects as well as the optical field tunneling process. Increasing dc bias F 0 would reduce the transmission peaks. One may notice that the emission current density increases even when W < ħω, whereas in Fig. 6c, the emission current density decreases when ħω /W > 1. The discrepancy arises because the relative value of the Fermi level E F /W (thus the normalized electron initial energy) changes when W changes in Fig. 7, but remains the same when the laser frequency changes in Fig. 6c.
It is clear that our exact analytical solution is derived based on the simple potential landscape defined in eq. (1). It is important to know the sensitivity of such calculations to the shape of the surface potential. The effects of more realistic surface potential have been studied previously, by solving directly the time-dependent Shrodinger Equation (TDSE) numerically, for example in refs 1 and 23. In Fig. 8, we compare our analytical calculation using potential landscape eq. (1) with the numerical simulations in ref. 23, which is for the more realistic potential is the image charge potential, the last term represents a finite Gaussian laser pulse of pulse length τ, and the remaining parameters are defined in eq. (1). In Fig. 8, the following parameters are used 23 , work function W = 5.5 eV, Fermi energy E F = 4.5 eV, laser wavelength is 800 nm, τ = 30 fs (~11 laser cycles), and dc electric field F 0 = 0.2 V/nm. For the potential profile without Φ m , our analytical results for the emission current density w vs. laser field F 1 are in excellent agreement with the numerical simulations. The inset in Fig. 8 shows the multiphoton order n as a function of Keldysh parameter γ, by fitting the slope of ∝ w F n 1 2 , again, showing very good agreement with the numerical simulation. These results serve yet another validation of our exact solution. More importantly, it indicates that our exact solution using infinite plane wave laser excitation is even applicable to more realistic laser profiles, as long as the pulse length is larger than the wavelength. With the inclusion of Φ m , the slope of the multiphoton regime is reduced in general, due to the reduction of effective potential barrier near the metal surface. It is easy to estimate that the potential maximum is when considering the static potential barrier only (i.e. no laser field), where W eff is the effective work function, which is found to be 4.96 eV for the given parameters. Similar approach for the effective work function has been adopted to study above-threshold photoemission 37 . By adopting this effective work function in our analytical calculation, we found fairly good agreement with the numerical simulation up to F 1 ~ 8 V/nm. Thus, by replacing W with W eff , our solution may still be used to estimate the electron emission characteristics in the multiphoton regime, even for potential profiles with image charge.

Summary
In summary, we have presented an exact solution to laser induced strong-field electron emission from a metal surface under dc bias. We analytically solve the time-dependent Schrödinger equation with both a single frequency oscillating field F 1 cosωt and a static electric field F 0 . The solution takes the general form of the superposition of electron plane waves with ladder-state energies (channels) separated by the photon energy ω  coupled to the electron's initial energy. It is found that increasing the dc bias F 0 would increase the electron emission, by opening up more channels, including the processes such as multiphoton absorption, optical tunneling, multiphoton emission, and field emission. When F 0 is small, clear transition from the multiphoton absorption regime to the optical field emission regime is shown. When F 0 is large, multiple processes contribute to the electron emission. In the limit of F 0 → 0, the solution approaches that of Yalunin et al. 20 , for the special case of F 0 = 0. In the limit of large F 0 , the solution recovers the Fowler-Nordhiem dc field emission. The time-dependent emission current reveals that intense current modulation may be possible even with a low intensity laser, by merely increasing the applied dc bias. The effects of laser frequency and metal work function are also examined in detail. It is found that for given values of electric fields, F 0 and F 1 , the maximum electron emission current density always occurs around the laser frequency ħω/W = 1. The emission current decreases with the laser frequency when ħω/W > 1. All electron emission mechanisms are captured in a single formulation, demonstrating transition from n-photon process to optical and field tunneling.
Our theoretical results exhibit good coincidence with previous experimental measurements on: measured kinetic energy spectra of the photoelectrons, for the transition from multiphoton regime to the strong field regime 17 ; shifts on peak of energy spectra 28,29,37 ; the scaling of voltage-and laser power-dependent measured electron flux 2 ; and periodic variation (at the laser period) of the electron yield due to the carrier-envelope phase (CEP) effects 17 . Thus, our exact analytical solution provides an efficient tool to analyze the complicated nonlinear dynamics of laser induced electron emission. It offers plausible interpretation of experimental results, and guidance of future experimental design. While our results are derived for continuous-wave excitation, they are shown to provide a good approximation to time-varying or pulsed fields (Fig. 8), as long as the optical frequency is much larger than the bandwidth of the driving field 4,20 . A slowly varying envelope approximation, then allows a straightforward extension of our results to pulsed excitation, as stressed in ref. 20.
In this formulation, we have made the following assumptions: 1) the electrons are emitted from a single energy level at the Fermi level; 2) the surfaces of the electrodes are flat and the problem is assumed one-dimensional; 3) the image charge potential is not directly included to make the analytical treatment possible. The effects of electrode geometry 40,41 , space charge 6,36,42,43 , nature of the ion lattice of the electrodes, thermal redistribution, and short pulse illumination may be considered in future studies.