Previous studies of terahertz-field enhancement have opened up many promising applications, such as optical switching devices [1] and highly sensitive molecular and biological detectors [2]. In particular, the huge field enhancement of up to 25,000-fold [3, 4] is truly remarkable, suggesting that an extremely high electric field could accumulate in an acutely small nanogap, of size as small as 1 nm. More importantly, recent studies have revealed that such field enhancement induces a strong field tunneling of electron wave packets [5, 6] in the vicinity of the strong electric field supplied by a hugely enhanced THz field inside the gap.

Even though a few works have discussed the relation between the incident THz field and the tunneling current associated with such strong field tunneling, a systematic theoretical study has been missing, though it is important for understanding subsequent experimental work on THz-field-induced electron tunneling. In this report, we show how those driving parameters such as carrier-envelope phase (CEP) difference of the incident THz field, strength of the electric field inside the gap due to THz-field irradiation, and external dc bias across the gap, affect the tunneling current measured across the nanometer-sized gap. We find that the phase is a key parameter in determining not only the total tunneling current, but also the manner of increasing the tunneling current with the other parameters, such as the incident THz field strength and the applied dc bias across the gap. We also suggest a novel method to determine the CEP difference, by measuring the tunneling current across the gap.

To understand and estimate how many electrons can be emitted in a nanogap with enhanced THz-field irradiation, it is worthwhile to introduce a simple picture in which we assume that the tunneling probability of the electrons at a given instant in time is simply dependent on the electric field strength associated with the incident THz field E(t), while following the Fowler-Nordheim tunneling formalism. Additionally, the incident THz field can be regarded as a dc field, because its time scale is 1000 times as long as the dwell time of electrons inside the nanogap. When the strength of the incident THz field is 10 kV/cm and the field enhancement factor is 1000, which is relatively easy to achieve using state-of-the-art technology in high-intensity THz generation and in extremely narrow nanogap fabrication [3], the electric field inside the nanogap reaches as high as 1 V/nm. If the gap size is about 5 nm, the total energy of an electron U_k obtained from this field is calculated to be as high as 5 eV. Here the terminal velocity of an electron accelerated by this field can be obtained as

From this calculation, the time of arrival (TOA) of the electrons at the opposite side of the gap should be less than 5 fs. Noting that a half period is ~500 fs for a single-cycled pulse of 1 THz frequency, our estimation is justified, and moreover the recoil or quiver motion of the tunneled electron is expected to be completely suppressed.

For more precise calculation of the electron tunneling current in various situations, a simulation method is developed by modifying previous models [7-10]. Here a time-varying THz field is assumed simply as

where f is a field enhancement factor and E₀, τ, ω, and ϕ are the peak electric field strength, pulse width, frequency, and carrier-envelope phase difference of the incident THz pulse respectively. Unlike in previous models, spatial variation of the electric field inside the gap is not taken into account, because of the extremely small gap size compared to the wavelength, which is ~λ/100,000 [3]. Those values are assumed to be τ = 1 ps and ω = 0.5 THz for every simulation performed in this manuscript. The emission probability p(t_B) for a given birth time t_B is obtained by applying the Fowler-Nordheim tunneling equation under the quasi-dc condition discussed above, i.e.

Here E(t_B) is the enhanced electric field induced by the incident THz field, Φ is the work function of the metal forming the walls of the nanogap, e is the elementary charge, m is the electron’s rest mass, and ℏ is the reduced Planck constant. By integrating p(t_B) over the duration of the incident pulse, the tunneling current can be obtained. When integrating, it is most important to note that electron emission can occur from either side of the gap. As schematically depicted in Fig. 1, for positive phase of the incident THz field, electrons tunnel from the left side of the gap (see Fig. 1(a)), but for the opposite phase, electrons tunnel from other side of the gap (see Fig. 1(b)). Hence the tunneling current for each phase destructively contribute to the total current, leading to a different tunneling current for the corresponding CEP of the incident THz pulse, which should be carefully considered in experiments. A theoretical consideration is supplied in the following section.

To calculate the terminal kinetic-energy distribution when an electron hits the other side of the gap, a semiclassical simpleman model is applied [7, 11]. In our situation it is important to know the exact TOA for each electron, so we observe the trajectories for those electrons emitted at each birth time, and find the time of impact t_n when the displacement of an electron becomes larger than the gap width of 5 nm. In this instance, the velocity of the electron v(t_B) is obtained as a function of birth time, and finally the terminal kinetic energy is found simply by applying . Rough estimation confirms that our electrons are in the nonrelativistic regime. For those electrons that move too slowly to reach the other side of the gap during the time window applied in the simulation, the terminal kinetic energies are taken at the instant when the simulation finishes; this treatment is valid when the simulation window is larger than the pulse width, because there would be no further acceleration after the THz pulse had passed.

Once the terminal kinetic energy is obtained, it is possible to calculate the kinetic-energy spectrum of the electrons when they hit the other side of the gap. The emission probability P(t_B) obtained above can be rewritten as follows:

where f(t) is a probability density function. What we want to have is a probability distribution P(E, E+ΔE) in the energy interval ΔE that can be given in the form

where f(E) is a probability density as a function of energy. Assuming that there is one-to-one correspondence between time and energy, i.e. E = g(t), the change of variable-gives the following relation:

with

Hence we have

Finally, the penetration depth of electrons should be considered while they travel inside the nanogap, because the typical nanogap is fabricated by inserting an ultrathin dielectric material, such as Al₂O₃ [12] or any of various two-dimensional materials [13]. In this simulation we assume that the gap is filled with Al₂O₃, of which the penetration depth x of an electron in microns is given as x = 0.1E^1.5/ρ [14], where E is kinetic energy of the incident electron and ρ is the density of the material, given as 3.89 g/cm³ for Al₂O₃. Considering that the width of the gap is 5 nm in our simulation, the corrected electron tunneling probability N_pen to successfully reach the other side of the gap can be given as

Hence, the tunneling current as a function of time J(t) has the form

and

As discussed previously, the effect of CEP on the tunneling current is primarily considered in our calculation. Figure 2 depicts the calculated kinetic-energy distribution as a function of birth time for different CEPs of 0 and pi/2 (Figs. 2(a) and 2(b) respectively), and the integrated tunneling current as a function of CEP through one pulse (see Fig. 2(c) for a linear plot and Fig. 2(d) for a polar plot). Here the incident pulse is assumed to have a duration of 1 ps and a frequency of 0.5 THz. The peak electric field is assumed to be 10 kV/cm, providing a field enhancement of 1000 (in the results, the enhanced peak electric field is 10 MV/cm, or 1 V/nm), which corresponds to the minimum electric field strength to generate the tunneling electrons in the simulation. As depicted in Figs. 2(a) and 2(b), the kinetic-energy distribution as a function of emission time almost resembles the electric field strength of the incident THz pulse (see the black curves in Figs. 2(a) and 2(b), having a maximum of ~5 eV for CEP of 0 and ~4 eV for CEP of π/2). This confirms our assumption discussed in the previous section. More importantly, the emission current J_T as a function of time is completely different for CEP of 0 and π/2, such that only a single peak near 0 ps is seen when CEP is zero, while two distinct peaks are seen when CEP is π/2. One should note that the polarity of the electric field is opposite for each phase, and thus the direction of electron emission is opposite, canceling out each other, and hence the net current J should be zero. Such behavior is well displayed in the plots of total tunneling current as a function of CEP (Figs. 2(c) and 2(d)). Here the tunneling current is almost zero near CEP of π/2 and 3π/2, a minimum at CEP of 0, and a maximum at CEP of π (see Fig. 2(c)). A polar plot (see Fig. 2(d)) shows clear dipole-like features (here the absolute value of tunneling current is used), which indicates that CEP is the most critical factor in tunneling-current generation.

Figure 3 depicts the kinetic-energy spectra for CEP of 0, π/2, π, and 3π/2, measured at the side opposite the emission site. It should be noted that the kinetic-energy spectra are highly monochromatic, as found in recent experimental work [15], regarding the slow variation of the THz field, which was discussed also in the previous section. Relatively high velocity allows electrons to reach other side of the gap before the carrier phase is inverted; hence, almost every electron presents only subcycle motion [11]. Additionally, the peak in kinetic energy lies at ~5 eV for CEP of 0 and π, and ~4 eV for CEP of π and 3π/2, suggesting that the acceleration of electrons can be controlled by CEP variation, which has been demonstrated in the near-infrared (NIR) region in previous work [16, 17]. Those results emphasize that the CEP of the incident THz field is critical in tunneling-current measurement, such that a detectable signal shows up only when the CEP is 0 or π; in contrast, no detectable signal is expected for π/2 and 3π/2. Those who try to implement this phenomena in an experiment should find a way to control, or at least to confirm, the CEP of their THz pulses.

In Fig. 4 we plot the total tunneling current as a function of amplitude of the incident THz field, varying the electric field strength from 16.7 kV/cm to 137 kV/cm with a fixed CEP of zero. As expected, the plot shows typical Fowler-Nordheim behavior of the kinetic-energy distribution and emission probability, for incident electric fields of 33.3 kV/cm (see Fig. 5(a)) and 137 kV/cm (Fig. 5(b)). Interestingly, the overall shapes of the kinetic energy distributions resemble each other, while the tunneling current as a function of time J(t) show different shapes. Notably, the full width at half maximum (FWHM) increases considerably, from ~100 fs to ~300 fs. Such broadening is attributed to the decrease in nonlinearity of the Fowler-Nordheim-type tunneling mechanism. The broadening also affects the kinetic-energy spectrum, as depicted in Figs. 5(c) and 5(d), such that the spectrum becomes less monochromatic at higher field strength, even though the maximum energy increases by a factor of ~8 (see Fig. 5(d)). These results suggest that for narrow emission one needs to decrease the incident THz pulse to be as small as possible, in both the time and frequency domains.

The observations and discussions above indicate that control of field emission of electrons in a nanogap by a strong THz field is a considerably delicate matter. Specifically, careful CEP control of the THz field (which conventionally is not considered) is required, in addition to careful intensity control, if one needs to obtain a highly monochromatic electron wave packet. To manage such delicacy, it is desirable to apply an additional dc bias, giving an additional electric field other than the THz field, for more efficient electron generation. The more important contribution of dc bias is a symmetry breaking in the nanogap, such that the electric field across the gap gives a positive effect only in a certain direction, where E_tot = E_THz +E_DC. Even a THz pulse having a CEP of π/2 could show different peak amplitudes for positive and negative phases, as experimentally manifested in previous work [10].

The manner of increasing the tunneling current as a function of the applied bias should also be different, with respect to the CEP differences. As schematically illustrated in Fig. 6, the amplitude of electric field with main emission phase where the peak of THz field is seen develops critically when an external bias is applied, for a CEP of 0, yet for a CEP of π the trend is completely opposite, such that the absolute value of peak electric field decreases with increasing dc bias. As a result, the bias dependence of tunneling current is quite different for different CEP values, as depicted in Fig. 6(c). Here the incident THz field and field enhancement factor are set at 10 kV/cm and 1000 respectively. As seen in the inset, the increase in tunneling current is much smaller for CEP values other than zero, being almost invisible using a linear scale for CEP of π/2, π, and 3π/2. The log-log plot shows the qualitative differences more clearly. In case of zero CEP difference, nonlinearity continuously increases with increasing bias. For CEP differences of π and 3π/2 (those plots completely overlap each other, due to symmetry), still such a trend persists, but the magnitude of the tunneling current is ~10⁴ to ~10² times smaller, compared to the CEP difference of zero. In contrast, for a CEP difference of π the trend is completely different, such that nonlinearity slightly increases at low bias of ~2 V and shows almost no change at other bias values, and even no detectable signal at bias values below 1 V. The result matches well with our schematic explanation that increasing dc bias destructively contributes to tunneling-current generation in this phase. All these results confirm the novelty and advantage of our method for determining the CEP difference of an incident THz field by observing the tunneling current as a function of dc bias, and examining the nature of the variation in nonlinearity.

To investigate further the combined effect of dc bias and THz field on tunneling-current generation, we plot the kinetic-energy distribution, tunneling current as a function time J(t), and kinetic-energy spectrum. Figures 7(a) through (c) depict the kinetic-energy distribution and J(t) measured at the dc bias of 2 V (0.4 V/nm), for CEP differences of zero (Fig. 7(a)), π/2 ((b), and π (c). The case of a CEP difference of 3π/2 is omitted, due to symmetry. It is clearly seen that the peak kinetic energy is much more evolved for a CEP difference of zero, compared to a difference of π. As a result, the peak in the kinetic-energy spectrum decreases for π/2 and π, compared to a CEP difference of zero (see Figs. 7(d)~7(f)). Specifically, a grey area where negative current is generated still persists for CEP differences of π/2 and π. In addition, J(t) shows a different behavior, in that peak values are 2 to even 9 orders of magnitude lower for CEP differences of π/2 and π, compared to a CEP difference of zero (see red curves in Figs. 7(a)~7(c)).

In conclusion, we have developed a theoretical model describing the tunneling behavior of electrons in an extremely small nanogap a few nanometers wide, when a THz pulse is incident, by modifying a conventional semiclassical model that is widely applied for near-infrared wavelengths. Using this model, we have investigated the effects of CEP difference and strength of the incident THz field on tunneling current across the nanogap, and found that the CEP difference is a critical parameter in net current generation along with the symmetry of the nanogap and the phase variation of the incident THz pulse. In addition, the increase in tunneling current follows the conventional Fowler-Nordheim equation. More importantly, we found that a dc bias also contributes to tunneling-current generation, but the character of this contribution is completely different for different CEP values. We expect that our results will be a good starting point and useful for further theoretical and experimental studies of strong-field tunneling phenomena in the THz frequency range.