Quasiclassical theory of non-adiabatic tunneling in nanocontacts induced by phase-controlled ultrashort light pulses

We theoretically investigate tunneling through free-space or dielectric nanogaps between metallic nanocontacts driven by ultrashort ultrabroadband light pulses. For this purpose we develop a time-dependent quasiclassical theory being especially suitable to describe the tunneling process in the non-adiabatic regime, when this process can be significantly influenced by the photon absorption as the electron moves in the classically forbidden region. Firstly, the case of driving by an ideal half-cycle pulse is studied. For different distances between the contacts, we analyze the main solutions having the form of a quasiclassical wave packet of the tunneling electron and an evanescent wave of the electron density. For each of these solutions the resulting tunneling probability is determined with the exponential accuracy inherent to the method. We identify a crossover between two tunneling regimes corresponding to both solutions in dependence on the field strength and intercontact distance that can be observed in the corresponding behaviour of the tunneling probability. Secondly, considering realistic temporal profiles of few-femtosecond pulses, we demonstrate that the preferred direction of the electron transport through the nanogap can be controlled by changing the carrier-envelope phase of the pulse, in agreement with recent experimental findings and numerical simulations. We find analytical expressions for the tunneling probability, determining the resulting charge transfer in dependence on the pulse parameters. Further, we determine temporal shifts of the outgoing electron trajectories with respect to the peaks of the laser field in dependence on the pulse phase and illustrate when the non-adiabatical character of the tunneling process is particularly important.

Abstract. We theoretically investigate tunneling through free-space or dielectric nanogaps between metallic nanocontacts driven by ultrashort ultrabroadband light pulses. For this purpose we develop a time-dependent quasiclassical theory being especially suitable to describe the tunneling process in the non-adiabatic regime, when this process can be significantly influenced by the photon absorption as the electron moves in the classically forbidden region. Firstly, the case of driving by an ideal half-cycle pulse is studied. For different distances between the contacts, we analyze the main solutions having the form of a quasiclassical wave packet of the tunneling electron and an evanescent wave of the electron density. For each of these solutions the resulting tunneling probability is determined with the exponential accuracy inherent to the method. We identify a crossover between two tunneling regimes corresponding to both solutions in dependence on the field strength and intercontact distance that can be observed in the corresponding behaviour of the tunneling probability. Secondly, considering realistic temporal profiles of few-femtosecond pulses, we demonstrate that the preferred direction of the electron transport through the nanogap can be controlled by changing the carrier-envelope phase of the pulse, in agreement with recent experimental findings and numerical simulations. We find analytical expressions for Together with quantum interference and entanglement, tunneling is one of the core phenomena characterizing the essence of quantum physics. For all these three phenomena comparison to the classical description benchmarks new possibilities opening in the quantum world. In the case of tunneling we have also a formalism connecting both the classical and quantum description, represented by the quasiclassical Wentzel-Kramers-Brillouin (WKB) method, first proposed in a general mathematical context of linear second order ordinary differential equations [1]. This formalism in many cases allows to obtain analytical or semi-analytical solutions and gain additional insights into their behaviour based on an extended classical intuition. The original quasiclassical approach is, however, suitable only in energy-conserving situations with static potential barriers or in some formally equivalent cases of time-dependent potentials which can be mapped to static descriptions since the time variable can be effectively seen as a spatial coordinate [2,3]. The seminal work of Leonid V. Keldysh [4] established a connection between the picture of tunneling and multiphoton ionization induced by laser fields, i.e. in temporally changing spatial potential barriers. This breakthrough achievement was followed by development of quasiclassical approaches applicable in situations of tunneling through spatial energy barriers varying with time [5][6][7]. Thus also the regime of so-called non-adiabatic tunneling [8][9][10][11], when a considerable energy is absorbed in the process of the underbarrier motion, could be captured within the same physical picture. Conceptually time-dependent quasiclassical approaches can be related to the path integral formalism [12], extended to the complex time plane and corresponding generally complex trajectories. In atomic physics, under certain conditions a formal derivation can be based on the strong-field approximation (SFA) [13][14][15][16]. Developments in ultrafast photonics leading to the appearance of tailored fewcycle laser sources opened new opportunities for studying of the control of pulseinduced dynamical tunneling from atomic systems [17,18], metal surfaces [19] and plasmonic nanoparticles [20]. It was demonstrated that the carrier-envelope phase (CEP) of the light pulses can be utilized as a control parameter that found its practical application, e.g., as a method of the CEP measurement [21][22][23]. It was further shown both experimentally and theoretically that strong ultrashort pulses can populate the conduction band in dielectrics [24] and semiconductors [25] due to nonresonant interband Landau-Zener tunneling of the electrons from the valence band into the conduction band and generate currents in unbiased systems on ultrafast time scales [26,27], whereby the current direction is controlled by the CEP of the pulse.
Recently we can also observe an increased interest to tunneling in nanosystems, since it can open new perspectives for ultrafast nanoscale devices. So the quantum tunneling regime was predicted and observed for plasmonic systems with nanoscale gaps and resulting plasmonic response properties were studied [28][29][30]. Another highly interesting case is realized in the tunneling microscope configuration where the charge transport between the tip and the surface can be strongly influenced by an external light field [31]. Moreover, application of tailored femtosecond few-cycle pulses (FCPs) in this configuration opened possibilities for sub-cycle coherent manipulation of the charge transport through the nanogap [32][33][34][35], enabling angstrom-scale spatial and subfemtosecond temporal resolution in tunneling microscopy [36]. In contrast to strongly spatially asymmetric configuration of the tunneling microscope, metallic nanoantennas with nanoscale gaps represent devices having a symmetric stationary part of the spatial potential barrier so that the whole control of the charge transfer across the gap can be realized solely via the CEP-controlled FCPs or their sequences [37,38].
In this work we develop a quasiclassical description of the charge transfer between two nanocontacts induced by CEP-controlled femtosecond light pulses and determined by the probabilities of the field-driven non-adiabatic tunneling. Therefore we consider the tunneling process in the corresponding potentials being both space-and timedependent and leading to values of the resulting Keldysh parameter γ, which determines the transition from the direct tunneling to the multiphoton regime [4], that can vary in a wide range. We base our theory on a Lagrangian formulation [7,[39][40][41] of the quasiclassical imaginary time method (ITM) [5,6] that brings, in our opinion, certain advantages with respect to the orginial ITM formulation concerning the justification of the derivation steps and finding the spatio-temporal structure of the resulting solution. With our approach we can address also the case of driving by the CEP-controlled FCPs that until now remained out of reach of the existing ITM considerations, where only solutions for even FCPs could be found [42][43][44][45], whereas a related approach of Keldysh provided us recently with a solution just for one case of a fixed-CEP odd single-cycle Figure 1. (a) Geometry of the nanogap (width d) formed by two metallic nanocontacts and the incoming driving few-cycle femtosecond light pulse with the electric field E(t) polarized along the z-direction. The nanocontacts are attached to an external electric circuit which allows to measure the electric charge transferred upon the pulse application. (b) The corresponding energy diagram for the electron travelling between the nanocontacts through the field-influenced time-dependent barrier. ∆E is the energy barrier for an electron at the Fermi level E F inside the metal. The red line shows an example how the energy of a tunneling electron changes with its position.
pulse [46]. Furthermore, the existing quasiclassical descriptions are rather suitable for free-space ionization problems whereas, as we will see, in the nanogap configuration the aspect of a small intercontact distance has to be addressed appropriately that eventually affects the solution structure.

General formalism
We will limit our consideration to one-dimension motion of the electron along the coordinate z. In the quasiclassical approximation the electronic wave function Ψ(z, t) is given by whereS(z, t) is the action. In order to find it as a function of z and t, one should find the general solution of the Hamilton-Jacobi equations [ The general solutionS(z, t) can be constructed using a complete solution where C is an arbitrary constant and L(z ,ż , t ) is the Langrange function, as Here C(t 0 ) is generally an arbitrary function of t 0 , whereas t 0 has to be expressed as a function of z and t using the equation Taking into account ∂S 0 The form of the function C(t 0 ) has to be selected in a way to match the boundary condition on the incoming part of the wave function at the position where the electron enters the region under the barrier, which we assume to be located at z = 0. At z = 0 and t = t 0 , the wave function is given by Ψ(0, t 0 ) ∝ e i S (0,t=t 0 ) = e i C(t 0 ) . Comparison of this expression with the wave function in the region outside the barrier in the neighborhood of z = 0, up to the exponential factors, allows to determine the form of the function C(t 0 ). If outside the barrier we can consider the electron energy E to be insignificantly influenced by the perturbation of the system, then we have and equation (6) turns into E = H t=t 0 . This is the case for an electron tunneling between two metallic contacts (see figure 1), because the electric field of the pulse is screened inside the contacts, due to a quick plasmonic response. Note that it is also the case for a charge carrier which is initially confined by a short-range potential [7, p. 295-297], [39,40]. The Lagrange function in equation (3) corresponding to one-dimensional motion of the electron in the laser electric field E(t) and a potential U (z) is given by where F (t) = eE(t) with e denoting the electron charge and m denoting its mass. z (t ) satisfies the equation of motion: with the boundary conditions Additionally, in this case from equation (6) we have where v 0 ≡ż (t = t 0 ) is the initial velocity and U 0 ≡ U (z = 0 + ) with z = 0 + being the position just beyond the barrier boundary [left and right limits of U (z) are different at z = 0 in the case of an abrupt potential step, which will be considered below]. In the case of the tunneling from the state with energy E this gives As far as for a tunneling state under the barrier ∆E < 0 [cf. figure 1(b)], we see that v 0 has a purely imaginary value. Solving the equation of motion together with the boundary conditions allows to determineS and v 0 =ż (t = t 0 ) as functions of z, t, and t 0 . Inserting the result for v 0 into equation (12), we can find then t 0 as a function of z and t. In result, the actionS can be expressed as a function of merely z and t. Thus, using equation (1), we are able to find the wave function Ψ(z, t). As we will see below, t 0 (z, t) is generally a multivalued function resulting in multiple solutions for Ψ(z, t). Within the presented formalism each of these solutions has to be analyzed separately and describes a part of the single tunneling process. However, there is no known recipe how to combine them together consistently. This issue is mostly relevant for the underbarrier dynamics close to the tunnel exit. Otherwise, the time-dependent probability density at each particular point in time t and space z is typically dominated by just one of these solutions (see section 3). If the electric field has several oscillation cycles, there is an additional multiplication of the number of the solutions, which are separated in time approximately by the oscillation period.
Let us now assume that the electron exits the classically forbidden region before it reaches the opposite metallic contact at z = d, excluding for a moment the case of very small spatial gaps between the contacts from the consideration. In order to determine the tunneling probability with exponential accuracy it is sufficient to know the value of |Ψ(z, t)| 2 at the set of points in the region behind the barrier where it has its maximum. Taking into account equation (1), we can see that these positions are determined by the condition ∂ ∂z ImS = Imp = 0 and therefore Equation (13) fixes the optimal classical complex trajectories [16] of the tunneling electron z opt (t). The value of t 0 must satisfy z opt (t 0 ) = 0. We see that for the optimal trajectories the electron velocity becomes real in the classically allowed region, which is intuitively expected. For other possible trajectories the velocity generally remains complex. From equation (2) follows ∂ImS(z, t) ∂t where we have a conditional partial derivative on the left hand side. This equation means that the imaginary part of the action does not change with time along any chosen optimal trajectory when the electron moves in the classically allowed region behind the barrier. Moreover, in order to calculate the corresponding imaginary part of the action, technically we can select any arbitrary value of time t on the real axis, even one corresponding to a time moment before the ultrafast laser pulse has arrived. This is possible because there is always a purely real classical trajectory z re (t) coinciding with the optimal complex trajectory for sufficiently large times after the laser pulse application, t > t ex . This real trajectory z re (t) corresponds to a fictitious electron coming from the classically allowed region to the potential barrier and then reflected back. It can be obtained from z opt (t), starting at t > t ex and propagating back along the real time axis (cf. also [49,50]). The consideration above becomes, however, inapplicable if the width of the potential barrier becomes so small that the electron is not able to exit out of the barrier before reaching the opposite contact. For brevity, we will call this situation small-distance scenario. Equation (13) does not generally apply for this case since |Ψ(z, t)| 2 may continue to change inside the barrier with increase of z up to z = d. To determine the tunneling probability with exponential accuracy we have to analyze the shape of |Ψ(z, t)| 2 at z = d. Then in order to find the optimal complex classical trajectory we need to formulate an appropriate replacement for equations (13) and (14). The optimal trajectory should lead to the maximum |Ψ(z, t)| 2 at z = d with respect to different possible values of time t so that in place of equation (14) we get ∂ImS(z, t) ∂t Denoting the particular time moment when this condition is satisfied as t E , in the case of the optimal complex trajectory we searching for the second of the boundary conditions (10) takes the form With exponential accuracy, the tunneling probability is given by whereS is evaluated at any final point (t, z) belonging to the trajectory z re (t), unless the small-distance scenario is realized when we have to take (t = t E , z = d).

Action and the propagating wave packet
Solution of the equation of motion (9) with the boundary conditions (10) gives where the functions V(t ) and Z(t , t 0 ) are defined as and the complex parameter v ≡ż (t = 0) is given by Now t 0 in equations (18), (21) and (22) is a yet unknown complex parameter that should be expressed as a function of z and t with help of equation (12). For that a system of two real equations Rev(z, t, t 0 ) + ReV(t 0 ) = 0 , Imv(z, t, t 0 ) + ImV(t 0 ) = ± 2 m ∆E (24) for the real τ 0 ≡ Ret 0 and imaginary part τ e ≡ Imt 0 of has to be solved. The selection of a positive or a negative sign in equation (24) leads to the same physical results. Without loss of generality, we may restrict our consideration to the positive sign. Explicit solutions can be found when the laser pulse shape F (t) is fixed.
Having t 0 (z, t), the actionS(z, t) can be determined using equations (3),(4),(7) and (8). If the the Lagrangian entering as the integrand in equation (3) does not have any singularity points as function of the complex time we may use integration by parts and writẽ Its imaginary part is then given by Inserting t 0 (z, t) into equation (26) and equation (27) we can determine the corresponding wave function. In particular, if it accommodates the real trajectory z re (t), we obtain the shape of the propagating wave packet. The probability density is given by |Ψ(z, t)| 2 ∝ e − 2 ImS(z,t) .
In order to discuss a more involved situation when the Lagrangian does have singularity points, let us define the standard path δ s in the complex time plane. This path starts at t = t 0 , goes firstly parallel to the imaginary time axis from the point (τ 0 , τ e ) to (τ 0 , 0) and then parallel to the real time axis to (t, 0), thus connecting t 0 and t, whereby not crossing any singularity point. For such a path equations (26) and (27) can still be used. Moreover, these equations stay valid for any other path δ such that the closed pathδ = δ s −δ does not encircle any singularity points. All these paths lead to the same result and are therefore physically equivalent. The case of a general choice of the integration path is analyzed in Appendix A for one particular example of the pulse shape F (t), which is also studied in section 3.1. For clarity, the consideration there is limited to optimal classical complex trajectories. We found that also for trajectories encircling singularity points equation (27) holds. Considering the closed pathδ mentioned above, it is helpful to introduce the winding number n j for each singularity point j and the total winding number N = j n j . One can show that the standard path δ s , having N = 0, leads to the maximal possible value of the probability, whereas any paths with N = 0 give smaller values. Physically, for a pulse-driven tunneling process the latter paths correspond to multiple reflections in the induced dynamic potential. The quasiclassical approach considered here is applicable only when the difference in probabilities between δ s and paths with N = 0 is very large (for more precise formulation, cf. Appendix A). Below we restrict our consideration to this applicability region, where it is sufficient to take into account only trajectories with N = 0. Outside of this region the utilized quasiclassical description is not readily justified and can lead to inconsistencies [cf. discussion after equation (64)].

Optimal complex trajectory and tunneling probability
As indicated above, in order to determine the tunneling probability with exponential accuracy it is not necessary to possess information about the whole wave packet of the electron. The knowledge of ImS(z, t) along the optimal trajectory z opt (t) (or trajectories, if there are several of them) in the region behind the barrier suffices. Therefore the task reduces to finding z opt (t), together with the related real trajectory z re (t), followed by calculation of ImS(z, t) for any real t and z = z re (t).
In this case the solution of the equation of motion, equations (18) and (19), should satisfy the following conditions: (i) Im[z (t = 0)] = 0, (ii)ż (t = t 0 ) = i 2 m ∆E, which are also valid for any complex classical trajectory, and additionally (iii) Imv = Im[ż (t = 0)] = 0, due to equation (13). For (ii) we have selected the plus sign in front of the square root, in accordance with the convention we decided to use for equation (24).
The choice of t = 0 for (i) and (iii) is convenient but not compulsory as far as any real value of t can be taken here. From (iii), we see directly that v is real. Then, taking into account ∆E > 0, these three conditions result in − vτ e + ImZ(0, τ 0 + iτ e ) = 0, Generally, these three equations allow us to find the three unknown real quantities v, τ 0 = Ret 0 , and τ e = Imt 0 . When they are determined we can express the imaginary part of the action through a real integral as With equation (17) this gives the tunneling probability. However, generally we may find multiple results for it if there are several physically different solutions of (28)-(30) for {v, τ 0 , τ e }, e.g., as illustrated in section 3.2. Let us again consider separately the small-distance scenario where the electron does not manage to exit out of the barrier into the classically allowed region before reaching the opposite contact. Proceeding as above but now taking into account equations (15) and (16), we obtain where condition (iii) originates from equation (15). It has a clear physical meaning that at the moment the electron reaches the opposite contact it has a real value of the kinetic energy. Note that this value has actually to be negative since we have assumed that the electron does not leave the classically forbidden region. The listed conditions result finally in a system of five equations, − τ e Rev + (t E − τ 0 )Imv + ImZ(t E , τ 0 + iτ e ) = 0 , Rev + ReV(τ 0 + iτ e ) = 0 , for five real quantities: t E , τ 0 , τ e , Rev, Imv. With these quantities determined, the imaginary part of the action at z = d and t = t E is expressed as replacing equation (31) for the considered scenario.

Ideal half-cycle pulse
In order to demonstrate in detail how the approach works for a particular case allowing for analytical results, let us consider the following pulse shape (cf. Refs. [6,42,43]): Here the parameter Γ, having the frequency dimension, determines the pulse duration. One can notice that the temporal integral over such a field does not vanish, and that is generally forbidden for light pulses propagating in the far field zone. However, one might, e.g., assume that the corresponding experimentally realized pulses possess merely weak oscillating or decaying tails [51][52][53]. The regions of the opposite polarity in respect to the main half-cycle assure that the integral of the field over the whole time axis converges to zero. However, the dependence of the tunneling probability on the electric field strength is highly non-linear. As a consequence, the impact of the weak tails of the pulse on charge transfer processes governed by tunneling is negligible. Hence models like equation (38) may be used. They can be considered for a qualitative understanding what happens during a single half cycle of a FCP. Inserting equation (38) into equations (20) and (21) we obtain Here we restrict our consideration to one branch of the multi-valued complex logarithm function: its principal value (for a more general situation, see Appendix A). Setting t = t in these equations we use them in equations (22)- (24). This leads to an equation for complex t 0 = τ 0 + iτ e (or a system of two equations for real τ 0 and τ e ): providing us with t 0 (z, t).
Here is a characteristic length and is a generalized Keldysh parameter [4,16,44,54] for the case of ideal half-cycle pulses.
As mentioned in the previous section, we get several solutions for t 0 from equation (41) for each pair of values of z and t (even though we have limited the consideration to the principal value of the complex logarithm function). The behaviour of these solutions in dependence on the position in the (z, t) plane is analyzed in Appendix B.
For the considered type of the driving field there are two physical solutions which give dominant contributions to the resulting probability distribution. Other solutions may be neglected. There is a branch point of order 1 located at (z = z c , t = 0). By going around this point in the plane (z, t) one of the two main solutions for t 0 transforms into another. Behaviour of the multi-valued functions τ 0 (z, t) and τ e (z, t), limited to the two main solutions, along the line t = 0 is illustrated in figure 2(a). Their topological structure close to the point (z = z c , t = 0) is similar to the Riemann surfaces for the real and imaginary parts of complex function √ ζ in the neighborhood of the branch point ζ = 0. We can separate the two main solution branches from each other in order to obtain distinct single-valued functions. Dependencies τ 0 (z, t) and τ e (z, t) are illustrated in figures 3(a) and 3(b) for one of the main solution branches (let us name it first solution here) and the region t > 0. The position and time are normalized by z 0 and Γ, respectively. The dependencies τ 0 (z, t) and τ e (z, t) for the other main solution branch (second solution) are shown in figures 4(a) and 4(b). As we will discuss below, the second solution dominates for the spatial region where the electron is under the barrier as well as in a certain vicinity after it may exit from this region. This solution, however, does not possess the corresponding real classical trajectory z re (t) and decays at larger  (53) and e − 2 ImS ≡ P e ≈ 5.7 × 10 −8 is given by equation (58). Black contour line in (b) corresponds to τ e = 0.752 < 0.7545 and encircles a region with even lower values of τ e than for the optimal complex trajectory. distances at all times t , for which then the first solution overtakes the leading role.
Using equations (22), (39) and (40) in equation (27), we obtain where the dimensionless function ξ(z, t) is given by with Everywhere , as it is determined by equation (41). The dynamics of the probability distribution following from equation (44)  Concerning the first solution, notice that the maximum probability for any fixed time moment t remains the same. It is also preserved along the optimal trajectory for real z and t. Thus if we view the temporal evolution of the spatial probability distribution shown in figure 3(c) as the dynamics of the wave packet of the emitted electron we should take into account that this wave packet is not normalized in a way that the particle number is conserved. Clearly, such a normalization is beyond the exponential accuracy of the applied method. It may be partly recovered going beyond the zeroth order of the saddle-point approximation used in the derivation of the standard imaginary time method based on the strong-field approximation [15,16]. Alternatively, the wave packet can be normalized to match the total probability resulting from its form at a fixed time moment, e.g. at t = 0. Looking at figures 3(a)-(c) we can see that along the real classical trajectory we have τ 0 (z, t) ≡ 0, τ e (z, t) ≡ const and ImS(z, t) ≡ const. ImS(z, t) reaches its maximum value on this trajectory. τ e (z, t) has there a conditional minimum value, under the condition that this value stays possible for any t, including t → ∞. In fact, on the real classical trajectory τ e = min z τ e (z, t) t→∞ . Smaller values of τ e (z, t) are possible for finite t -see the region encircled by the solid black line in figure  3(b), but they disappear as t grows. In contrast to the first solution, the probability distribution corresponding to the second solution has an evanescent character away from the barrier, as can be observed in figure 4(c). The values of the probability at each time moment increase along with the electric field of the applied pulse so that the probability density is instantaneously dragged in the direction of the opposite contact.
To find the corresponding optimal complex trajectory, assuming that the electron exits into the classically allowed region before reaching the opposite contact, we use equations (39) and (40) for t = 0 and write the real and imaginary parts at this time moment explicitly: Let us find solutions with the absolute value of Γτ e being below π/2 (for other possible solutions as well for a clarification when and why they can be neglected, see Appendix A). Inserting equations (49) and (50) into equations (28)- (30) and eliminating the variable v leads to two equations for the determination of two real quantities τ 0 and τ e : Generally, such a problem should be treated numerically. Here the solution is simplified by the fact that the left hand side of equation (51) vanishes only if either τ 0 = 0 or τ e = 0. For τ e = 0 equation (52) cannot be satisfied. Thus we must have τ 0 = 0 and therefore f (Γτ 0 , Γτ e ) = cos 2 Γτ e . Then from equations (29) and (49) follows v = 0, i.e. the electron leaves the classically forbidden region with zero velocity. Notice that this property is not assumed but does follow from the derivation. As a result equation (52) simplifies to tan(Γτ e ) = γ HCP .
For fixed values of γ HCP and Γ this equation has only one solution in the range Γτ e ∈ (−π/2, π/2) giving τ e = 1 Γ (arctan γ HCP ). Having obtained τ 0 , τ e and v, we can then determine the full optimal complex trajectory: In case of the standard path δ s in the complex time plane, we can divide the trajectory into two parts: (1) underbarrier motion t where changes from iτ e to 0 and (2) the following motion in the classically allowed region where t is real and t > τ 0 = 0. For part (2) we find that z opt > w holds, where is the distance travelled by the electron under the barrier. Considering equation (54) only for real t but extending it to negative values, we get the real trajectory z re (t) of the fictitious electron reflected from the barrier (see figure 3). It is interesting to notice that the value of w, somewhat counterintuitively, does not coincide with the distance to the branch point z c . For example, for the parameters of figures 2-4 we have w = 0.3166z 0 whereas z c ≈ 0.3461z 0 , i.e. the difference is around 10%. This means that the second solution still dominates over the first solution not only within the dynamic tunnel barrier but also in its immediate neighbourhood. However, the probability density of the second solution is evanescent at larger distances and therefore does not induce there a density flow towards the opposite contact, in contrast to the first solution. These findings are in agreement with the numerical observations for the electronic density in the gap between contacts obtained by the time-dependent density functional theory (TDDFT) [38] in the case of relatively large gaps (∼ 6 nm). Thus our theory is able to contribute to the understanding of the pulse-induced tunneling dynamics of the electron density in vicinity of the nanocontacts. Furthermore, we will discuss below that in the strongly non-adiabatic regime, with large values of the Keldysh parameter, z c can considerably exceed w [cf. inset to figure 5(a)] increasing the role of the second solution, in particular for gaps 1 nm. In terms of finding optimal complex trajectories and corresponding probabilities, this situation is captured by the small-distance scenario, treated separately in the end of this section.
For the imaginary part of the action acquired along the optimal trajectory we obtain from equation (31) with v = 0 and τ 0 = 0: With equation (17) and the solution of equation (53) for τ e , we get the tunneling probability This is the result obtained in [42][43][44] by the ITM, where the problem of atomic ionization by an ultrashort light pulse was considered. In contrast to these works, here we use an alternative formulation of the quasiclassical approach to tunneling in time-dependent fields [7,39,40]. With that, firstly, all the intermediate steps have been consistently explained avoiding ad hoc assumptions. As we will see below, this is especially important for the following consideration of more realistic waveforms of the applied light pulses. Secondly, the details of the two main solution branches for the time-dependent, real-space quasiclassical wave function have been clarified, with important physical consequences highlighted below. Let us discuss the applicability range and limit cases. The quasiclassical description used here is appropriate only if condition ω ∆E is fulfilled for the frequencies ω belonging to the spectral content of the pulse. This implies here Another applicability condition is given by a general requirement that the tunneling probability remains small: (a) (b) Figure 5. (a) Function G(γ) entering condition (61). Inset shows the dependence of the distance travelled by the electron under the barrier w given by equation (56) and critical distance z c when γ HCP is varied in the same way and for the same conditions as described below. (b) Dependence of the tunneling probability on the Keldysh parameter. Green line shows the case of an ideal half-cycle pulse, equation (58). We start from the parameters as in figures 2 and 3, assume γ = γ HCP and then vary γ by altering the value of the peak electric field strength. Gray vertical bar with label "v" indicates the validity threshold (65). To the left of this bar we cannot rely on the result represented by the green line. Blue line illustrates the case of a few-cycle pulse with a Gaussian envelope with the temporal shape given by cos(ωt)e −t 2 /(2σ 2 ) , whereas we selected σ = 1/ω = Γ, connecting the time scales to the case of an ideal half-cycle pulse. Red line shows the case of continuous wave (CW) driving, also with 1/ω = Γ. In both latter cases, the calculation was based on the results given in [44]. Full black line corresponds to the approximation given by equation (62)  Finally, in our derivation of equation (58) we have assumed that the distance travelled by the electron under the barrier w does not exceed the distance between the contacts d. This restricts possible values of the generalized Keldysh parameter from above by the condition Failure to fulfill this condition always means that the small-distance scenario is realized, with a different procedure to find the optimal complex trajectory, as we described in section 2.3. We will return to this scenario below. The dependence G(γ HCP ) is shown in figure 5(a). For a fixed value of γ HCP , it follows that d should exceed G(γ HCP ) ∆E/F 0 . For example, we have G(γ HCP ) = 1/2 if γ HCP ≈ 1.585. This implies d > ∆E/(2F 0 ): During the motion inside the classically forbidden region the electron absorbs an energy being equal to ∆E/2. We will revisit conditions (59)-(61) in discussion section, deliberating on realistic system parameters.
In the direct tunneling limit, i.e. γ HCP → 0, we have w = ∆E/F 0 . No energy is absorbed in the process. Inequality (61) is satisfied when the barrier height ∆E does not exceed the potential drop across the intercontact region corresponding to the peak value of the electric field of the light pulse, i.e. d > ∆E/F 0 . Then decomposing the left hand side of equation (53) in Taylor series we get Γτ e ≈ γ HCP and therefore τ e → 0 with γ HCP → 0. In the leading order in γ HCP in the exponent, equation (58) simplifies to which is just the well-known expression for the quasiclassical probability of tunneling through a static triangular barrier. From equation (62) and equation (60) follows a condition limiting the electric field strength ‡: For γ HCP 1, the applicability of the quasiclassical approach would break down for smaller values of F 0 than dictated by equation (63).
In the multiphoton limit, with γ HCP → ∞, formally we get Γτ e → π/2 and equation (58) reduces to where ∆E/( π −1 Γ) is the average number of absorbed photons in the multiphoton transition. Notice that in contrast to the conventional continuous wave case here photons of various energies belonging to the pulse spectrum participate in the process. However, it is immediately clear that equation (64) cannot represent the correct limit case result because this limit can be reached by decreasing the electric field amplitude and keeping other system parameters constant. Obviously, when the electric field vanishes the probability must also vanish that is, however, not the case for equation (64). In fact, as we discuss in Appendix A, the utilized quasiclassical description becomes invalid for too large values the (generalized) Keldysh parameter. The validity is restricted by the condition In order to obtain the full picture, it is important to calculate the optimal complex trajectories and corresponding probabilities also in the small-distance scenario. Looking at equations (32)-(32), we can eliminate there variables Rev and Imv by expressing from equation (34) and ‡ In the case of atomic ionization this corresponds to the requirement that the maximum applied electric field should be still much smaller than the characteristic atomic electric field. from equation (35). Then we can find from equations (39) and (40) − i arg (cos Γτ e cosh Γτ 0 + i sin Γτ e sinh Γτ 0 ) that by inserting into equations (32), (33) and (36) leaves us with a system of three equations for three variables: t E , τ 0 , τ e . Solving this system we can determine these variables and find then the corresponding probability. However, we can notice that equation (36) actually means that Imv = 0 or/and Rev + ReV(t E ) = 0 must be fulfilled. The first option would lead actually again to equations (51) and (52) and then to equation (53), obtained above for the case when the electron leaves the classically forbidden region before reaching the opposite contact, that would contradict to the assumption of the small-distance scenario. Choosing the option Rev + ReV(t E ) = 0 and using equations (68) and (69), equations (32), (33) and (36) can be recast as In general the system of these three equations has to be solved numerically but we can notice that there is always a solution with τ 0 = 0, leading also to t E = 0. This results in Γτ e [γ HCP − tan(Γτ e )] − ln cos(Γτ e ) = d z 0 , from which we can find then also τ e . Note that this equation actually coincides with equation (41) if we take there z = d and t = 0. Thus we have actually already analyzed its solutions and illustrated them in figure 2 for a particular choice of parameters. For d < z c there are two solution branches, whereas it is actually the second solution which plays the dominating role in this scenario. Moreover, for w < d < z c where the first solution possesses also the optimal complex trajectory exiting into the classically allowed region before reaching the opposite contact and merging there with the real classical trajectory z re (t), the second solution, having the optimal complex trajectory without this property, still leads to higher values of probability. Other possible solutions of equations (70)-(72), which can be found numerically and cannot be related to our first and second solution branches, lead to much smaller probabilities than the dominating solution branch and may be therefore neglected. Finally, using equation (37) we can find so that the resulting tunneling probability in the small-distance scenario can be expressed as Dependence of the tunneling probability on the Keldysh parameter is illustrated in figure 5(b) for the studied half-cycle pulsed driving in comparison with the driving by a FCP with a Gaussian envelope, continuous wave (CW) driving as well as the result that follows from equation (62) (direct tunneling limit). The validity threshold following from equation (65) is indicated by the gray bar labeled by "v". Note that for large electricfield strengths, and correspondingly small values of γ, the results for the tunneling probability converge to the direct tunneling limit in all cases. On the double logarithmic scale chosen for figure 5(b) it leads to an exponential behaviour log 10 P = −ae ln10 log 10 γ (a > 0). This agrees with the corresponding limit case result given in [14] (cf. p. 1797 there) but seems to disagree with some later numerical calculations [31,36] used to model related experimental data and showing linear behaviour on the double logarithmic scale, log 10 P = −a − b log 10 γ (a , b > 0), in the limit γ → 0. However, the latter calculations are actually based on the Reiss theory [14] so it remains unclear to us why they fail to reproduce the result of the direct tunneling limit inherent to the Reiss theory (or generally Keldysh-Faisal-Reiss theory) as mentioned above. The experimental data provided in [31] would not allow to differentiate strictly between the two dependencies. Concerning the experimental and theoretical data in [36] we do not recognize the underlying physical reason for a quite abrupt slope change at γ ∼ 1 indicated there and cannot confirm it by our theory.
In discussion of the behaviour of the probability for large values of γ HCP the crucial aspect is to know the value of the intercontact distance d. For very large d the probability is always determined by the first solution and decays rapidly with increase of the Keldysh parameter, though somewhat slower than for the CW driving and even more slower than given by the direct-tunneling expression. In this consideration we have to keep in mind the restriction imposed by the condition (65). However, due to the decrease of the probability to very small values the discussion of the first solution for extremely large values of γ HCP has anyway no practical sense. Moreover, for d 1 nm it is the second solution that plays the dominating role with much larger probabilities for large γ HCP , not the first solution. At finite d the second solution has a clear physical limit case for vanishing field strengths, i.e. for γ HCP → ∞ in the context of figure 5(b). This limit case corresponds to tunneling via a static rectangular potential barrier of height ∆E and width d, as we also illustrate in figure 5(b) for d = 0.6 nm. We see that at such small gaps the second solution overtakes the first one already at values γ HCP below 1, determined by the condition z c (γ HCP ) = d, and leads then to much higher probabilities than the first solution for larger values of the Keldysh parameter. In result, the whole dependence of the tunneling probability on the field strength, converted to the Keldysh parameter, is significantly different from the result taking into account the first solution only or equivalently the assumption of large d. This has to be kept in mind when discussing related experiment results similar to whose presented in [31] and [36].

Realistic few-cycle pulse
To model a realistic FCP [51] we use a wave packet with a rectangular spectral shape [55] and a variable flat phase. The temporal profile determined by where the Fourier transform F (ω) is constant and can be decomposed into a real positive amplitude F 0 and a phase factor factor e iφ : Here ∆ω = ω 2 − ω 1 is the spectral bandwidth of the pulse and φ is the CEP. In this part we limit our consideration to the optimal complex classical trajectory and tunneling probability, keeping in mind the structure of the quasiclassical wave function discussed in the preceding subsection. Using equations (76) and (77) in equations (20) and (21) we get Substituting these expressions into equations (28)- (30) and eliminating the variable v leads to the following equation system for the times τ 0 and τ e : with and i n (x) denoting the modified spherical Bessel functions of the first kind [56]. For each pair τ 0 and τ e we find the optimal complex trajectory and the tunneling probability determined by where The behaviour of the solutions of equations (80) and (81) in dependence on the CEP is illustrated in figures 6(a) and 6(b) [for the temporal shape of the chosen FCP, see figure 7(a)]. We see multiple solutions in these figures. Each solution can be actually attributed to a peak of the FCP. Tunneling takes place at a time moment close to the temporal position of the corresponding peak, in the sense that there is a value of τ 0 being close to that position (but generally not exactly coinciding with it, as we will discuss below). The absolute value of the imaginary tunneling time τ e increases monotonically with the distance of the peak from the center of the pulse as it is changing with φ. In figure 6(c) for each of the solutions we show the tunneling distance w travelled by the electron under the barrier, whereas for the case when this quantity is smaller than the distance between the contacts. We see that w behaves non-monotonically in dependence on the phase. For 6 main solutions, having the smallest values of |τ 0 | and |τ e | with respect to the remaining solutions, the values of w are confined to a relatively narrow interval, which for the selected parameters constitutes around 0.05 nm.
For a FCP the distance between contacts d with respect to w brings an additional aspect for the resulting charge transfer. When d is large enough with respect to w the electron might have no chance to reach the opposite contact after the tunneling through the barrier because the electric field changes its polarity after a certain time and might drive the electron back to the original contact. The corresponding motion is determined by the corresponding classical trajectory z re (t). In figure 7(a), as an example, we show two such trajectories originating from two positive peaks being closest to the center of the pulse at φ = π/2. Whereas the trajectory belonging to the right peak does not return back to the left contact, the electron appearing in the gap around the left peak moves a certain distance towards the right contact but then is turned back due to the changed polarity of the field and flies back to the left contact. If the distance between the contacts is too large the latter electron never reaches the opposite contact. Then, neglecting possible reflections by the contact boundaries, it does not contribute to the resulting charge transfer. For larger gaps between the nanocontacts it is thus essential to model not only the tunneling process but also the following dynamics in the classically allowed region. An example can be found, e.g., in [38]. Here we focus our attention on the regime when d > w for any of the contributing trajectories but d is still small enough so that the tunneling process determines the charge transfer. Before that let us look more closely at the positions of τ 0 associated with the exit point of the standard complex classical trajectory to the classically allowed region.
In figure 7(b) we zoom into the time interval close to the position of one of the field peaks. We can observe that the corresponding value of τ 0 does not exactly coincide with the position of the peak. In the figure we also show the electron velocity found from the classical trajectory z re (t). The time moment where the velocity vanishes determines the turning point of the trajectory. We see that this time moment generally does not coincide neither with τ 0 nor with the peak-field position. This fact is remarkable in view of recent suggestions to associate the tunnel exit time with the time moment when the velocity vanishes that can be determined, e.g., by the backpropagation of the emitted electronic wave packet [49,50,57]. One justification argument behind these suggestions is that "the pure tunnelling dynamics in a semi-classical description is characterized by an imaginary momentum component in at least one degree of freedom" whereas "all momenta are real in classically allowed regions". However, in general the ad hoc assumption that the momentum of the electron during its underbarrier dynamics has a vanishing real part, implied by the ITM, has no justified reason to hold. As it is in fact happens for the situation of figure 7(b), this momentum runs generally over complex and not purely imaginary values. Thus at the moment when the electron enters the classically allowed region and its momentum becomes real, it should not necessarily also vanish. In fact, it depends on the particular chosen path in the complex time plane. For the standard vertical path, when t changes from t 0 to Re(t 0 ) = τ 0 , considering the situation of figure 7(b) we can see that the momentum does not vanish at t = τ 0 but acquires a finite real value. The momentum, of course, vanishes if an alternative, non-vertical path reaching the real time axis at t withż re (t) = 0 is selected. The discussion on the definition and meaning of the tunneling time with the related controversy [54,57,58, and references therein] is out of the scope of the present paper but we hope that our findings bring an additional useful insight in this context.
It is interesting to analyze the attosecond time shift ∆τ 0 between τ 0 and the temporal position of the corresponding peak of the electric field. We illustrate this in 7(c). For φ = 0 there is no shift for the central negative [in terms of the induced force F (t) = eE(t)] peak, whereas the electron is emitted later (earlier) when it might be expected for the proceeding (succeeding) peak. When the CEP φ is changed away from zero the time shift ∆τ 0 appears also for the central peak, due to the break of symmetry. Generally, the absolute value of ∆τ 0 grows as the position of the peak moves away from the center of the FCP, with its sign being positive (negative) for τ 0 < 0 (τ 0 > 0). It would be interesting to find if this effect is possible to measure experimentally.
Finally, having obtained the solutions for t 0 based on equation (83) we can calculate the resulting transition probabilities belonging to each of these solutions. In order to represent the result for all main solutions simultaneously it is convenient to fix a particular peak of the electric field and the related solution and then evaluate the dependence of the tunneling probability on the phase in the extended phase scheme, where the CEP φ can take any real values, not limited by the interval (−π, π). The result is shown in figure 8(a) for two different values of the peak electric-field amplitude. For comparison we also plot the results corresponding to the quasi-static approximation when the probability is calculated using the direct-tunneling approach with the static potential determined by the value of the electric field at the corresponding peak. We see that the quasi-static approximation significantly underestimate the probability, especially for lower field values. For lower F 0 the effect is pronounced already at φ = 0 (central peak). The non-adiabatic enhancement grows with increase of the absolute value of the CEP, as the distance between the peak and the FCP centre raises. Therefore, in order to determine the charge transfer from the whole FCP it is sufficient to take into account just several contributions coming from the closest peaks to the FCP centre, with the respective sign determined by the direction of the field.
In order to calculate the resulting probability difference between positive and negative contributions we limited the consideration to the six closes peaks to the FCP centre, with τ 0 and τ e limited to the intervals indicated in figures 7(a) and 7(b). For each electron emerging from the tunneling region we calculate its classical trajectory and take its contribution into account only if it reaches the opposite contact avoiding in the meantime the original contact. For larger values of the distance between the contacts d it happens that the classical field-induced dynamics plays a major role in determining the overall charge transfer [38]. It is then insufficient to limit the consideration to optimal complex trajectories z opt (t) and related real trajectories z re (t). All complex trajectories (or in other words whole emitted electronic wave packets) emerging at each field peak should be taken into account. It can happen that the part of the wave packet reaching the opposite contact does not contain the optimal trajectory at all. In order to have a situation dominated by the tunneling process we consider a configuration when d is close to the tunneling distance w for all real trajectories but still always exceeds it by a certain amount, so that we may consider only the wave-packet-like (first) solution (cf. section 3.1) neglecting the evanescent (second) solution. In this regime we can also neglect the effect of interference between electronic wave packets emitted at neighbouring electricfield peaks [16,44]. We take the parameters of figure 6(c) and select d = 0.75 nm.
The resulting probability difference is shown in figure 8(b). For comparison we show also the results for the quasi-static approximation when the electron would not gain energy during tunneling. As expected, we see that the non-adiabatic description leads to significantly higher values of the probability difference than the naïve quasistatic approach. We can also nicely observe the modulation of our calculated quantity determining the total charge transfer induced by the FCP in dependence on the CEP. Thus the direction of the charge transfer can be controlled by FCP on an ultrashort time scale, in agreement with the experimental observations and other theory predictions [36][37][38]59]. In the studied case the modulation can be fitted by a cosine function so that the difference of the fit to the calculation result in figure 8(b) is visually barely distinguishable and we therefore did not plot it separately.

Discussion of the relevant system parameters and approximations
Let us at first briefly discuss the relevant parameters of light pulses (cf., e.g, reference [51]). Taking typical pulses in the near-infrared with wavelength λ ∼ 1.1 µm we have the half-cycle duration ∼ 2 fs that gives the relevant scale for 1/Γ in case of the considered ideal HCPs. Realistic FCPs can be as short as τ FWHM ∼ 6 fs. The peak electric-field strength E 0 can be estimated from the pulse energy E p ∼ 200 pJ that is focussed to an area of A ∼ 20 µm 2 . Taking into account that the pulse energy can be expressed as E p ≈ c 0 2 |E 0 | 2 Aτ FWHM , where c denotes the speed of light and 0 is the vacuum permittivity, we arrive at E 0 ∼ 1 V/nm. Between the nanocontacts the field is amplified due to the plasmonic enhancement by a factor θ, which can take values up to ∼ 100 but it is strongly dependent on a particular realized configuration. The pulse typically also experiences certain phase shifts or/and distortion. In our modelling we operate with the anticipated field in the gap. In our illustrating examples, we orient ourselves on E 0 = 4 V/nm and E 0 = 10 V/nm inside the gap.
Next, we want to review the parameters of the nanocontacts. We assume that they are made of gold having the Fermi energy E F ∼ 5.5 eV. Depending on the manufacturing method, currently different sizes of the nanogap width d are possible: from around 30 nm [37] to 6 nm [38,59] and further to subnanometer values in break junctions [36,60]. The height of the effective energy barrier is determined by the gap medium or it can be also influenced by the properties of the utilized substrate. For example, we would have ∆E ≈ 5.1 eV for a Au/vacuum/Au composition of the nanocontacts and ∆E ≈ 4.2 eV for a Au/SiO 2 /Au junction. In our calculations we took ∆E ≈ 5 eV.
Based on the above parameters of the driving light and nanocontacts we can estimate other important relations for our study. The relation of the average photon energy to the energy barrier width amounts to Γ/∆E ∼ 0.07 so that one of the conditions for the validity of the quasiclassic approach given by equation (59) can be seen as well satisfied. Concerning the (generalized) Keldysh parameter, assuming E 0 = 4 V/nm we estimate γ HCP = 0.94. This value corresponds to the intermediate regime between the direct tunneling and multiphoton ionization and we used it for many of our illustrations. For E 0 = 10 V/nm we have a 2.5 times smaller value of γ HCP . The other validity condition given by equation (60) is also well satisfied for all parameter values discussed above.
Further, it is useful to review the characteristic spatial scales of the investigated system. The maximum tunneling distance for an electron at the Fermi level w F,max = ∆E/(|e|E 0 ) constitutes ∼ 1.25 nm taking E 0 = 4 V/nm. It decreases to just ∼ 0.31 nm for E 0 = 10 V/nm. For the electrons at the bottom of conduction band in the contacts the maximum underbarrier distance to overcome would be ∼ 2.6 nm (∼ 0.65 nm) for E 0 = 4 V/nm (E 0 = 10 V/nm). The characteristic length z 0 defined by equation (42) is estimated to ∼ 2.8 nm (∼ 0.7 nm) for E 0 = 10 V/nm (E 0 = 4 V/nm) with 1/Γ = 2 fs.
Beyond the evaluation of the tunneling probability for a single electron at the Fermi level we can roughly estimate the total number of electrons transferred between the nanocontacts by the applied pulse. We take the electron density of gold n Au = 5.9 × 10 28 m −3 , assume that the area of the nanocontacts constitutes ∼ 100 nm 2 , calculate the Fermi velocity of the electrons v F from E F as 1.4 × 10 6 m/s, and estimate the number of electrons with energy close to the Fermi level and hitting the boundary in the temporal vicinity of the field peak as (100 nm 2 ) × (2 fs/10) × (n Au /20) × v F ∼ 100. Depending on the resulting probability difference being in the range 10 −7 −10 −3 this gives 2 × 10 −5 − 2 × 10 −1 transferred electrons per pulse. With the standard pulse repetition rate of 40 MHz this amounts to ∼ 10 3 − 10 7 electrons per second or approximately 10 −4 − 1 pA. Note that whereas our theory does not allow for a quantitatively precise evaluation of the magnitude of the transferred charge it provides its dependencies on various parameters of the system, and that in an analytical or semi-analytical way.
There is one effect neglected in our consideration that potentially can significantly influence the tunneling barrier and resulting probability. When the electron is outside the metal in a static case it should experience the interaction with its own image charge created inside the metal contact. Taking the image charge into account leads to the modification of the triangular shape of the barrier leading to the so-called Schottky-Nordheim barrier [61,62] given by where the factor 16π appears here in place of the usual 4π for the usual Coulomb potential because the distance from the tunneling electron to the image charge is the double of the distance to the metal, leading to the additional factor 1/4 in the force and therefore also in the potential. This leads to the lowering of the maximum height of the potential barrier with respect to ∆E by [63] For E 0 = 4 V/nm (E 0 = 10 V/nm) this amounts to ∼ 2.4 eV (∼ 3.8 eV) reduction and could even remove the barrier completely for higher fields. However, we have to keep in mind that we deal here with a charge transfer process that takes place on an ultrashort time scale. In fact, the charge density reorganization leading to the appearance of the image charge interaction requires some time to be formed, determined by the inverse plasmon frequency. The image charge effect becomes essentially time-dependent and can be taken into account, e.g., by using the corresponding velocity-dependent potential [64,65]. One can expect that the resulting barrier reduction is considerably lower than following from equation (86) [64][65][66][67]. We think that our quasiclassical approach, due to its time-dependent nature, has a good potential to be able to incorporate the image charge effect and to ultimately clarify the importance of the dynamic barrier reduction.

Conclusion and outlook
We have presented a quasiclassical theory for the description of tunneling and charge transfer in nanocontacts that is driven by half-cycle and phase-controlled few-cycle pulses. The theory is capable to account for the dynamics of the underbarrier motion and the energy absorption taking place during this process. Based on a simple model of an ideal half-cycle we are able to construct analytical solutions for the main solution branches of the electronic wave function in the classically forbidden region as well as after exiting out of this region. We have derived the expression for the tunneling probability that was already known from the imaginary time method in the case of sufficiently large distances between the contacts but now it has been determined for any intercontact distances. We have found that for larger intercontact distances the solution branch corresponding to a rising electron density inside the barrier and forming an outgoing wave packet plays the dominating role in most of the classically allowed region whereas the solution branch with the falling electron density inside the barrier eventually forms an evanescent wave that can be neglected, unless we consider the behaviour in vicinity of the tunnel exit. Here the evanescent-wave solution already may start to give a larger probability. This is especially pronounced in the strongly non-adiabatic regime with higher values of the Keldysh parameter, where there is an extended spatial region where the evanescent-wave solution dominates over the solution corresponding to the outgoing wave packet with a classical trajectory. It is one important finding of this work that this effect has to be taken into account if the boundary of the second contact occurs to be in this region. At very short distances such that the electron stays always in the classically forbidden region the roles of the solution branches definitively interchange and the falling-density solution branch plays the leading role. Based on these results we were able to calculate the dependence of the tunneling probability on the strength of the applied field, also in the regime where a crossover between both solution branches takes place. Further, in case of the phase-controlled few-cycle pulses we used our theory to find analytical solutions for the complex tunneling times and probabilities which determine the amount of the total induced charge transfer through the nanogap. In particular, we studied the configuration when the intercontact distance is such that the outgoing wavepacket solution with corresponding real trajectories can be used to calculate the resulting probabilities. We compared the turning points of these trajectories with the real parts of the tunneling times and with the temporal positions of the field peaks. We found that in general all these quantities are different, with temporal shifts being in the attosecond range for typical parameters. The amount of the transferred charge and the transfer direction can be controlled by the carrier-envelope phase of the pulse, whereas the values obtained in the non-adiabatic regime are significantly higher comparing to the quasi-static direct tunneling approximation.
One natural extension of the studied problem can be a consideration of a more complex spatial structure of the tunneling region, e.g., having an additional quantum well inside the tunneling barrier [68]. In such a case one may expect an intriguing interplay between the tunneling process steered by an ultrashort pulse and the energy level structure inside the quantum well that should enable an ultrafast selective population of the levels. In order to treat such a problem, a generalization of our method beyond the exponential accuracy [16,69] is probably required, which also represents an important and interesting task by itself. Finally, an appropriate inclusion of the dynamical image effect [64][65][66][67] into the description might lead to further improvement of our understanding of the light-induced tunneling in nanocontacts.

Appendix A. Paths encircling singularities in the complex time plane
For the pulse shape as in equation (38), Z(t , t 0 ), defined in equation (21) and entering equations determining t 0 , is a multivalued function. The value depends on a particular integration path δ between t 0 and t . For real t , we can calculate Z(t , t 0 ) as a sum of an integral along the standard path δ s and an integral along the closed pathδ = δ s − δ, as discussed in the end of section 2.2. Proceeding in this way we write The integral in the second term can be evaluated using Cauchy's residue theorem and we get where the total winding number is a sum of all winding numbers N (δ, a k ) of the pathδ for each singularity point a k = π 2 + kπ (k ∈ N). The remaining integral along δ s in equation (A.2) has two parts: one wheret runs parallel to the imaginary axis from (Γτ 0 , Γτ e ) to (Γτ 0 , 0) and another from (Γτ 0 , 0) to (t , 0). This leads to Evaluating the second integral in the brackets we obtain In figure A1 we illustrate the behaviour of the solutions of this system in dependence on n tot for the case when n = 0. Moreover, for this figure we restrict Γτ e to Γτ e ∈ (−π/2, π) for a better overview of possible paths corresponding to these solutions.
In figure A2 solutions of equations (A.12),(A.13) are shown for a larger range of Γτ e , including also the possibility of n = 0. Selecting an arbitrary solution, we indicate a possible path corresponding to this solution so that its characteristic features are clear. For paths that may encircle singularities in the complex time plane, in place of equation (31) at first we have to use (A.14) However, at least for the pulse shape given by equation (38), we can find that the two terms in the curly brackets of the integral on the second line of equation (A.14) lead to contributions having opposite signs and eliminating each other. Therefore, equation Calculating the probabilities P for various possible t 0 in the complex time plane based on equation (A.15), we found that all values of t 0 with negative imaginary part τ e lead to P > 1. All these solutions have probability densities rising with the distance inside the barrier. They also have a property that the electron initially moves to the left, away from the barrier, when moving in the complex time plane from the corresponding t 0 directly towards the real axis. Such solution branches are unphysical in the context of the posed tunneling problem and therefore can be ruled out from our present consideration.
For solutions in the upper part of the plane, i.e. with τ e > 0, we found that the maximum value of the tunneling probability corresponds to the solution with n tot = n = N = 0 having the standard path δ s as a possible integration path in the complex time plane (see blue path in figure A1). Typically with increase of the values of |n tot | or n the probability values drop extremely rapidly with respect to the standard-path solution. In such a situation we can justify neglecting all solutions with non-vanishing n tot or/and n . The situation changes if the respective probabilities do not vary drastically with respect to the standard-path solution. This occurs, e.g., when the amplitude of the driving electric field is decreased, resulting also in the correspondingly increased value of γ HCP . In figure A3 we illustrate how the probabilities corresponding to solutions from the part of the complex time plane shown in figure A1 behave in dependence on n tot when the amplitude of the driving field is varied. We can see that for larger electric fields, with γ HCP ∼ 1, the change of n tot from 0 to 1 already leads to a drop in the probability value as large as several hundred orders of magnitude. However, when the electric field is decreased and the value of the probability of the standardpath (n tot = n = N = 0) solution declines towards the value given by equation (64) the probabilities of n tot = 0 solutions in contrast rise. For small enough fields, which for our choice of parameters in figure A3 correspond to γ HCP 100, the dependence of P (n tot ) in the neighbourhood of n tot = 0 saturates to the constant level given by equation (64). In this limit there are many solutions delivering comparable values of the probability. However, within the utilized quasiclassical description there is no straightforward way to combine the corresponding multiple solution branches taking interference effects into account. Since we cannot use solely the standard-path solution in this case, the validity of the whole approach in its form presented in the current paper breaks down. Where would we then set the validity threshold in terms of the value of γ HCP ? Looking at figure A3 we can argue that this threshold is observed when the deviation of the standard-path probability given by equation (58) from the value of the limit-case probability given by equation (64) is comparable with the latter value. Using for the argument of the exponential in equation (64) that for large γ HCP we have arctan γ HCP = π 2 − arctan(1/γ HCP ) ≈ π 2 − (1/γ HCP ), we obtain the validity condition (65). Note that this subtle issue had remained insufficiently clarified in references [6,[42][43][44], whereas it was actually addressed in the last published paper of Leonid V. Keldysh [46] based on the Keldysh-Reiss-Faisal approach, comparable with other quasiclassical descriptions. Moreover the approach of [46] allowed to obtain appropriate limit case expressions for the tunneling probability. We agree with Leonid V. Keldysh that [46] "the weak field regime seems to be of more academic interest: for such short pulses, the effect is hardly experimentally observable." However, we should remark that for us it has been also important to resolve any apparent unexplained paradoxes, like finite probability values in the limit of vanishing strengths of the driving field following from equation (64), to ensure the overall consistency of the utilized method. Finally,   Figure A3. Tunneling probabilities for solutions with Γτ e ∈ (0, π) (cf. figure A1) in dependence on n not and the electric field amplitude scaled by a dimensionless factor E n . The probability is determined from equations (17) and (A.15) with γ HCP = 0.94/E n and F 2 0 /( mΓ 3 ) = 34.7E 2 n , correspondingly. The inset shows the magnified region of larger probabilities where log 10 P is closer to zero. The black horizontal line in the inset indicate the limit-case value for vanishing electric field after equation (60).
one should also mention that whereas the issue connected to multiple solution branches, originating physically from a multiple-reflection behaviour of the electron moving in the dynamical potential induced by the field, occurs for the pulse shape given by equation (38) it may be just absent for other pulse shapes [43,46]. are shown in figure B1(a). The behaviour of −|h(Γt 0 )| in the plane τ 0 , τ e at each of these points is then illustrated in the corresponding figures B1(b)-(f).
As we noticed in Appendix A, the solutions in the lower part of the complex time plane are unphysical and may be ignored. In figures B1(b)-(f) we observe that there are two main solutions, which are imaginary for small z. With increase of z these solutions move towards each other along the real time axis until they merge at z = z c . Then with further increase of z they again split and start to move parallel to the real time axis away from each other [see also figure 2(a)]. Having the same imaginary part of t 0 , these two solutions then lead to the same value of the probability at t = 0 and large z