Analysis of a higher-energy structure in nanotip enhanced fields

We investigate strong field ionization of an atomic gas in a plasmonically enhanced field resulting from the illumination of a nanometer-sized structure with ultrafast laser pulse. We use perturbation theory to derive an approximate solution for electron’s motion following ionization. These analytical estimates are corroborated by the time-dependent Schrödinger equation and classical trajectory Monte Carlo simulations. Notably, our approach can be used to obtain electron energy spectra without having to rely on numerical simulations. This allows for a deeper study of the dependence of electron energy spectrum on the properties of the near-field, suggesting electric field sensor applications. We derive an analytical expression for the location of the peak of the higher-energy structure (HES) as a function of laser parameters and near-field decay length. We find a particularly strong dependence of the energy peak on laser frequency, with lower frequencies causing a significant upward shift in the final electron energies. Combined with control of the width of the HES, which can be done by changing the size of the nanostructure, this points to the possibility of using nanotips as sources of ultrashort electron beams of tunable energy.


Introduction
Part of the 2018 Nobel Prize in Physics was awarded to high intensity ultrashort laser fields. The peak intensity of these fields is high enough to distort the binding energy of atoms and molecules, leading to a number of well-studied effects in strong field ionization and related phenomena [1]. For instance, ultrafast tomography relies on the diffraction of the ionized electron from the ionic core of the parent ion [2]. The most fundamental results involving strong field interactions rely on the above-threshold ionization (ATI) and the closely related phenomenon of high-harmonic generation (HHG) [3][4][5][6][7][8].
While high intensity fields, in the range of 10 14 W cm −2 , are necessary for the above described phenomena, it is possible to achieve them using only moderate strength oscillators plasmonically enhanced by nanotips, which could dramatically enhance the local electric field intensity by orders of magnitude relative to the driving field [9]. While relying on more accessible laser technology, this approach also introduces a wealth of new and important phenomena due to the spatial inhomogeneity introduced by the plasmonically enhanced fields [9][10][11][12][13][14][15].
When a nanotip is exposed to light with wavelength much longer than its radius, the surface electric field is significantly enhanced and the near field decays exponentially. Nano-optics, also called near-field optics, is based on this property and enables a range of applications such as in microscopy, spectroscopy and in ultrafast electron imaging [9,10,13,15,16]. The inhomogeneous near field plays an important role in the electron movement after ionization, significantly affecting the final electron energies. Recently HHG has been reported in crystalline targets by using locally enhanced fields [12,[17][18][19], and extensive theoretical works have been conducted to investigate the plasmon-enhanced [20][21][22][23][24] HHG. An extension of the harmonic cut-off was observed, which may lead to the production of extreme ultraviolet light and generation of attosecond pulses. At the same time, experimental and theoretical work demonstrates a significant extension in electron energy spectrum due to plasmonically enhanced fields [11,[25][26][27]. While most theoretical work has focused on the gas in the vicinity of a nanostructure, a similar analysis can be applied to solids.
In strong field ionization of atomic gas, a dipole approximation is typically used, which neglects the spatial dependence of the laser field. In this case, the ATI energy spectrum exhibits a cut-off around 2U p for direct emission, and a cut-off of 10U p from elastic back-scattering, where U p is the ponderomotive energy [3][4][5]. However, in the vicinity of nanotips, the dipole approximation typically breaks down due to considerable spatial variation introduced by the exponentially decaying near-field [28]. When this spatial variation is significant relative to the electron excursion, the electron trajectory is considerably altered [29], while theoretical modeling becomes considerably more complex. This photoemission regime, where a transition from quiver to sub-cycle motion takes place, is described by the parameter δ, which is determined by the ratio of the plasmonic field decay length to the electron excursion [30,31].
The spatial dependence of the inhomogeneous near-field can be derived by solving the Maxwell equations incorporating both the geometric and material properties of the nanotip and the applied laser pulse characteristics [32]. There are several numerical approaches for simulating the near-field, such as boundary element method [33][34][35], finite-difference time-domain and finite element simulation [36][37][38]. A semi-classical approach for obtaining the near-field distribution relies on an analytical multipole expansion method [39]. The numerically retrieved electric field can be approximated by a power series [21,27,32] and, in some cases, by a linear dependence [22,32] and a trapezoidal envelope [21].
Recently, an ATI photoelectron spectra from atoms in inhomogeneous field in the tunneling regime was numerically studied using classical trajectory Monte Carlo (CTMC) simulations combined with a solution of a time-dependent Schrödinger equation (TDSE) [40]. The electron energy spectrum exhibited a pronounced higher-energy structure (HES) in the region of 2U p . Using CTMC simulations, it was shown that the HES originated from the direct electrons ionized within a narrow time window of the laser pulse [40]. Overall, this showed [40,41] that the HES can be explained in terms of classical electron trajectories in the presence of the laser field and enhanced near-field. This provides a promising method to link ionization time with final energy, and study the influence of field inhomogeneity, laser intensity at the location of atoms, laser wavelength and the number of cycles inside the pulse on the HES.
In this work, we study the HES from ATI of atoms using perturbation theory applied to a two-step model, where the first step is ionization and the second step is the classical description of electron's motion. We obtain the HES distribution and derive the HES peak position, demonstrating that the energy distribution shifts from a monotone decrease in the homogeneous field to a more complicated behavior in the inhomogeneous field. We furthermore investigate the dependence of electron energy distribution on a wide variety of field parameters, including pulse length, field inhomogeneity, laser wavelength, and intensity. We discover additional features of the HES, characterized by secondary peaks near the main peak. Note that this is the first work that obtains final electron momenta distributions analytically, independently of TDSE or CTMC simulations.
The rest of the paper is organized as follows: in section 2, we briefly describe the linear approximation of the near field and the Ammosov-Delone-Krainov (ADK) ionization probability distribution. Then we derive an analytical expression for the HES peak position. In section 3, we study the dependence of HES distribution on different parameters, including field inhomogeneity, laser intensity, laser wavelength and the number of cycles in the laser pulse. We compare the HES features with analytical expression and explain the origin of the HES. We summarize our results and present conclusions in section 4. Throughout this paper, the intensity used in theoretical calculations means the intensity experienced by atoms in the near field, unless otherwise indicated.

The nanotip near field and classical description of electron motion
The electric field enhancement is strongest near the tip apex of the nanostructure, which can be described using a hemisphere of radius R. Away from the surface, the near field decays exponentially and relaxes to the strength of the incident field. The decay length, l f , scales approximately linearly with R. We consider an atomic gas in the vicinity of a nanotip illuminated by the laser light. Using a linear approximation, the near field distribution can be expressed as follows [42]: where cos 2 ωt 2N is the carrier envelope containing N laser cycles per pulse, φ 0 is the carrier envelope phase (set to φ 0 = π), andx is the laser field polarization direction. The inhomogeneity parameter is given by β, which is half of the inverse of decay length 2β = 1/l f . The atoms are placed outside the leftward-pointing nanotip, such that the electric field increases toward the nanotip. The sign of 1 + 2βx is not reversed for any electron trajectory considered.
The electron dynamics in the laser field can be obtained by solving the classical motion equation, It has been demonstrated that the solution of equation (2) and the equation of motion without Coulomb potential typically deviate little [41]. Hence, we neglect the Coulomb force in the electron propagation. The nontrivial motion driven by the laser field along the x direction can therefore be approximated as One can solve this equation numerically as in the CTMC method, or one can use perturbative methods to analyze it and express electron displacement as a sum of successively higher-order terms. The details of the perturbative expansion are described in section 2.3.
No mater which method is used, initial conditions are required to solve the above equation. In our two-step model, electrons are first ionized at time t 0 by tunneling (calculated by ADK formula as in section 2.2), secondly driven by the inhomogeneous laser field to the detector. For the classical motion (second step), an ionized electron is assumed to appear at tunneling exit with zero velocity at time t 0 , carrying a probability w(t 0 ) and accumulating energy while driven by the laser field. In this paper, we plot the energy spectra of ionized electrons just like a parametric equation, each final energy is paired to a probability by the ionization time t 0 . One will notice in the figures that some energy points have multiple-valued probabilities, which are linked to different ionization time t 0 .

The ionization probability distribution
The first step in the two-step model is ionization. In the tunneling regime, where γ 1, the ionization rate of the atom is [43,44], where The laser field F (t) is given by equation (1), and I p = 13.6 eV is the binding energy of the target hydrogen atom, γ = I p 2U p is the Keldysh parameter, and Q c (γ) is the Coulomb correction of the tunneling rate, given by where n * = Z 2I p is the effective principle quantum number of the atom, Z is the charge of the atomic core.
The incident laser pulse can not only ionize the atoms (after enhancement), but also has a chance to ionize the metal nanotip itself. One question immediately raised on how to distinguish electrons ionized from atoms and nanotip? In fact, it is not necessary to worry about it, since electrons available for ionization in atoms and nanotip are not of the same order. The key parameter β used to describe the inhomogeneity of enhanced field, is determined by the decay length, thus related to the size of nanotip [40]. For the value of β in this paper, the size of nanotip is on the order of 10 nm, which could contain 5000 atoms on its surface, but the typical number of gas atoms in laser focus is about 10 10 [15]. The huge number difference makes the ionization from nanotip itself negligible.

Perturbation approximation
The second step in the two-step model is propagation, in the process of which the dynamic parameters can be deduced by solving the classical equation of motion.
Since the inhomogeneity parameter, β, is small even for very sharp nanotips (for example, β = 0.002 for l f = 12.5 nm), we can use the following perturbative expansion for electron displacement, where x 1 is the first-order, x 2 the second order and x i the ith order approximation. To solve for each successive order, we plug equation (7) into equation (3). We can then describe each order of velocity and displacement by integrable functions: And where v i is the ith order velocity correction. x int is the electron tunneling exit position, representing an initial condition for classical propagation and given by [45,46], where Z eff is the effective nuclear charge that the active electron sees during ionization. The current discussion involves ionization of neutral atom, the asymptotic behavior of it would be Z eff = 1.
Using ionization time, we can establish a one-to-one correspondence between the electron final velocity (given by the electron velocity at the end of the laser pulse) and the ionization probability. This is of particular importance, since the HES distribution we discuss later is determined by the final electron velocity, which takes the form:

The HES structure
The electron's final velocity can be calculated using perturbation approximation. Including the first two terms, we write, v (1) the corresponding final energy is It is known that the HES electrons are ionized near the peak of the laser cycle [40]. For the CEP = π laser pulses, this corresponds to t 0 = 0. The time integration starts from t 0 = 0 to t = Nπ/ω, when the cos 2 pulse reaches zero. To lowest order, the final velocity and displacement can be expressed as From equation (16), v 0 = 0, since N is an integer. Therefore we use the next order term to describe the relation between each parameter and the HES peak position. The next order term in the perturbative expansion is given by, for the case of N = 2. Then the peak position of HES can be expressed as Hence we have derived an analytical expression, such as equation (20) for N = 2, which gives a relatively accurate and simple approximation for the peak position of HES. In the following section, we use this expression to demonstrate how the field inhomogeneity and laser parameters affect the location of HES.
The series expansion for the HES structure is convergent, and keeping velocity terms up to third order, resulting in HES energy given by, We can also estimate the HES width by calculating the full-width at half-maximum (FWHM) of the HES distribution. This can be obtained by the energy difference of electrons liberated at two times near t 0 = 0 which carry half maxima ionization probability. We should note that CEP = π and N = 2 are used when analyzing the peak position and FWHM of the HES. For other values of CEP and N, the calculation can be modified while still following the procedure outlined above.

Results and discussion
The break-down of the dipole approximation due to plasmonic near field enhancement dramatically alters the electron trajectory following ionization. Consequently, the energy distribution profile is significantly changed, resulting in the HES, as well as other unique features. The analytical expression we derived for the HES peak position includes the inhomogeneity parameter, β, as well as laser parameters, such as intensity, wavelength and the number of cycles. Below, we compare our analytic result, given by equation (20), with the peak position obtained using equation (22), corresponding to the maximum of ionization probability. We also compare our results with the previously obtained approximate solution in [41]. We end with a discussion of the physical origin of the HES electrons.

Inhomogeneity scaling and the minimum width
Using perturbation theory, we have derived the expression for final electron energy, given by equation (22), as detailed in section 2. In the upper panel of figure 1 we show a series of HES spectra for β between 0 to 0.003 and fixed laser parameters, obtained by solving equation (22). Note the excellent agreement between our analytical results and the CTMC simulations implemented in [41].
In the case of a homogeneous field, the electron count decreases monotonically with energy, with the lowest energy electrons coming from ionization at the peak of the laser field. By contrast, in plasmonically enhanced fields, the electrons ionized at the peak get accelerated to higher energies, making up the HES. Overall, the entire electron spectrum shifts to higher energies and is no longer monotonically decreasing, as shown in the top panel of figure 1.
In the bottom panel of figure 1, the analytical expression for the peak position given by equation (20) (the black curve) is compared with CTMC (the blue squares), and the position obtained from equation (22) (the dark red circles), corresponding to the maximum of ionization probability. The CTMC points (the blue squares) were obtained following [41]. The four peak positions, at corresponding values of β, show a surprisingly good agreement between the two analytical solutions obtained using equations (20) and (22), and the numerical solution using CTMC.
It is worthwhile to note that our result shows a scaling law as a function of inhomogeneity, laser parameters and the HES characteristics. This scaling law is also applicable for lower laser intensities. Figure 2 shows the HES curves and the peak positions with laser intensity I = 2 × 10 13 W cm −2 and wavelength λ = 4400 nm, in order to keep the Keldysh parameter the same as in figure 1. The HES peak positions clearly show a quadratic dependence on inhomogeneity parameter β, which is consistent with the scaling equation (20).
The black and pink points in the upper panel of figure 1 represent the energy of electrons coming from before (black) and after (pink) the peak of the laser pulse, when the probability of ionization is half of the maximum probability. These two points overlap in the homogeneous field because of the symmetric  (21). Laser intensity at atoms I = 1 × 10 14 W cm −2 , wavelength λ = 2000 nm and the pulse duration N = 2. The black rhombus points indicate the half maximum ionization probability from the laser rising edge of main peak, while the pink square points indicate that in the falling edge. Bottom: the peak position of HES (E Kin ) vs the inhomogeneity β in the two-step model produced by our peak position formula (the curve), given by equation (20). The blue squares show the peak position obtained by CTMC and the dark red circles show the position obtained by the peak position from the upper panel. The inhomogeneity parameter is given in atomic units and the energy is given in eV. distribution in the final velocity. However, for small β, the electron coming from the falling edge of the laser cycle (pink) gains more energy than the electron coming from the rising edge (black). Interestingly, this situation reverses itself for large β value, leading to an intermediate case where the energy of HES electrons is sharply focused, as shown by the green curve in figure 1. This change clearly happens between β = 0.002 to β = 0.003. Figure 3. FWHM of energy distribution as a function of β. This is obtained by the final energy of electron ionized with half maximum ionization probability from the rising edge minus that from the falling edge and this lasts around 0.67 fs. FWHM first decreases to negative values then increases, and crosses zero around β = 0.0021, which indicates the narrowest energy distribution. We chose the same laser parameters as in figure 1. The red dotted line is a quadratic fit, which demonstrates the quadratic dependence of FWHM on β. The inhomogeneity parameter is given in atomic units and the FWHM is given in eV.
The nearly mono-energetic electron accumulation, indicated by the green curve in figure 1, has important practical applications, as noted in [40]. To investigate this further, we calculated the FWHM of HES for a series of β ranging from 0 to 0.003. Figure 3 shows the FWHM as a function of β, obtained by taking the absolute value of the difference between the half-maximum value coming from the rising edge, corresponding to black dot in figure 1, minus the half-maximum value of the falling edge, corresponding to pink dot. The actual difference between the two half-maxima, which takes negative values for β < 0.0021, is shown in figure 3 alongside the always positive FWHM curve.
Note that the FWHM intersects zero for homogeneous (β = 0) field, because the final velocities have a symmetric distribution around the peak of the laser field, with electron density decreasing monotonically as a function of energy. This is well-known in the context of strong field ionization of atomic gas. Hence, the sharpest accumulation of nearly mono-energetic non-zero energy electrons occurs around β = 0.0021, where the FWHM approaches zero. The difference between the two half-maxima is well-described by a quadratic fit, which shows that FWHM has a quadratic dependence on β.

Laser intensity scaling and wavelength scaling
The excellent agreement between the analytical perturbative approach introduced here, and the CTMC simulations motivates further analysis of the peak position dependence on laser parameters. In particular, this analysis can be more comprehensive with analytic formulation, where the functional dependence on all parameters is known, as compared to direct numerical simulations of electron dynamics, such as TDSE or CTMC.
The upper panel of figure 4 shows the energy resolved photoelectron spectra corresponding to four laser intensities (see caption for details) ranging from I = 0.8 × 10 14 W cm −2 to I = 1.5 × 10 14 W cm −2 for β = 0.002, λ = 2000 nm and N = 2. The HES shifts to higher energies as the intensity increases, as can be seen from the upper panel of figure 4.
The lowest order displacement, x 0 , has a linear dependence on F 0 , as can be seen from equation (17). The integral of F 0 × x 0 therefore has quadratic dependence on the laser electric field. Hence the peak position has quadratic dependence on laser intensity, as can be seen in the lower panel of figure 4. Careful readers may notice the small difference between the curve and marked circles when the laser intensity is I = 1.5 × 10 14 W cm −2 . This is introduced by the x int term, which appears in x 0 and v 1 , and plays a more obvious role in stronger fields. Figure 5 shows the dependence of final electron energy on wavelength. The wavelength range is in the 1500 nm-3000 nm, corresponding to δ range from 4.3212 to 1.0803. Interestingly, just as for the β parameter, there is an intermediate value of λ for which the HES structure is particularly narrow.
The peak position has a highly nonlinear 6th power dependence on λ, consistent with equation (20). This dependence is confirmed in the lower panel of figure 5. Since a quiver amplitude is given by l q = F/ω 2 , laser wavelength is crucial in determining whether the system falls in the intermediate regime, determined by the ratio of the field decay length in the vicinity of the nanotip to the quiver amplitude: δ = l f /l q = l f ω 2 /F.

The sub-peaks of the HES distribution
Prior subsections investigated the influence of inhomogeneity, laser intensity and laser wavelength on the HES distribution for short pulses, corresponding to N = 2. Here we discuss longer pulses, where the other cycles contributed significantly to the overall electron counts, as a result of a smaller ionization probability  difference between the main peak and the sub-peaks of the laser pulse. For longer pulses, there is no one-to-one correspondence between ionization time and final energy, since electrons coming from different cycles can contribute to the HES structure. This is in contrast to short pulses where the HES electrons come from a single narrow time window around the peak of the central cycle.
The contribution of different laser cycles for longer pulses is quantified in the upper panel of figure 6, which shows the HES distribution in the case of N = 5. The lower panel of figure 6 shows the laser field as a function of time, with different colors corresponding to the main peak (purple) and the two nearby sub-peaks (blue and green, respectively). The contributions from the sub-peak immediately after the main Figure 6. Top: the respective contributions of the main peak and two nearest sub-peaks to the electron counts as a function of final energy. Bottom: the laser electric field, including a pulse envelope, color-coded to indicate the time-window where the main peak and the two sub-peaks plotted above originated. N = 5, laser intensity I = 1 × 10 14 W cm −2 , wavelength λ = 2000 nm and β = 0.001. peak (green) is strikingly different from the sub-peak preceding the main peak (blue) due to electric field inhomogeneity in the vicinity of the nanotip.

Conclusion
In conclusion, we investigated the dependence of the HES on laser and nanotip parameters, such as wavelength, field intensity, decay length and pulse duration. We focused mainly on the intermediate regime, δ ∼ 1, where the electron quiver amplitude is too small to immediately escape from the near field but high enough to feel its decay over a single period of oscillation [30]. As noted in prior studies [40,41], approximate solutions in this regime are particularly difficult, and one normally has to rely on numerical simulations, such as CTMC or TDSE. Here, we have shown how to treat this regime analytically within perturbation theory.
In particular, we derived an analytic formula for the location of the peak of the HES. This formula shows quadratic dependence both on the field inhomogeneity parameter, β, and the peak intensity, as well as strongly nonlinear (6th power) dependence on the laser wavelength. Moreover, for the first time to the best of our knowledge, we obtained the electron spectrum for the intermediate δ regime analytically, without having to resort to numerical simulations. Our analytical results show excellent agreement with previously published CTMC simulations.
Using analytically obtained HES curves, we found that the width of the HES structure depends strongly on the nanotip size (which determines the near-field decay length, l f ), in agreement with prior numerical studies [40,41], as well as the laser wavelength. This suggests that the width of the HES peak has a strong dependence on the nanotip 'non-adiabaticity' parameter, δ ∝ l f ω 2 . We plan to explore this further in future studies. The dependence of the width and peak location of the HES structure on laser parameters and field decay length is of particular importance in future application. These applications involve the use of the HES as a sensor of optical near field and for production of localized nearly mono-energetic ultra-short electron beams for time-resolved electron-diffraction experiments.