Coulomb blocking of sequential tunnel ionization in complex systems

When atoms and molecules are exposed to intense low frequency laser fields, the dominant response is sequential tunnel ionization of charge states with increasing ionization potential. Sequential ionization is assumed to proceed as separate one electron processes. The theoretical analysis developed here reveals that in complex systems sequential tunnel ionization can be inhibited by Coulomb blocking. When ionization potentials of subsequent charge states are close to each other, multiple tunneling events can occur during a half cycle and in close proximity, so that a tunneled electron can block the next tunneling electron. In sub-nm clusters driven by near infrared single-cycle pulses, Coulomb blocking reduces two-electron sequential tunneling by up to 2-3 orders of magnitude.

Strong laser field physics presents the foundation for various research areas; among those are attosecond science [1,2], laser driven generation of highly charged complex systems for influencing chemical reactions [3,4], laser machining of micro-and nano-scale structures [5][6][7], the generation of ultrashort electron pulses [8][9][10][11][12] for imaging [13], and light-driven electronic nanocircuits [10,12]. Commonly, the key strong field processes, such as optical field ionization [14], high harmonic generation [15], and above threshold ionization [16], are treated as single-electron processes. Electrons are ionized sequentially and do not interact with each other during ionization and subsequent propagation in the continuum. Two notable exceptions have been identified so far. The first multi-electron process is nonsequential tunnel ionization [17]. The tunneled electron quivers in the continuum and, on re-collision with the parent system, kicks out another bound electron. The second process is shake-up/shake-off ionization [18,19], where the tunneling electron shares energy with the remaining electron core and excites or ionizes another electron.
Our theoretical analysis reveals the importance of another strong field multi-electron process. Strong field physics has been focused to a large extent on exploring the dynamics of atoms and small molecules, for which the above picture applies. Recently, more complex systems, such as clusters and nanotips, have come into focus [3,10,12]. In these systems the ionization potentials of the various ionic states are spaced more closely than in atoms. This increases the likelihood of ionizing more than one electron during a half cycle. Due to its resultant proximity, a tunneled electron can notably increase the barrier for the next tunneling electron, thus Coulomb blocking (CB) its exit. While CB has been studied in the context of nano-structures and -devices, such as tunneling junctions [20,21], its role in optical field ionization has not yet been considered.
The quasi-classical, adiabatic approach developed here is limited to two electrons tunneling from a model cluster. In the case of near-IR wavelengths, we find that CB appears most pronounced in smaller clusters, driven by single-cycle laser pulses with moderate field strengths, in which case electrons have little opportunity to avoid each other. Suppression of sequential two electron ionization by 2-3 orders of magnitude is found. This is of relevance for a number of single-cycle laser processes, including laser driven generation of single attosecond electron pulses from nanosurfaces, laser generated attosecond electronic signals in nanocircuits, and laser formation of highly charged complex systems, where ultrashort pulses are used to minimize heating and fragmentation. Potential experiments to measure CB are discussed at the end.
We start by introducing the quasi-classical, quasi-static approach used to model CB of tunnel ionization. The complex system is represented by a spherical model cluster. Pauli blocking is not considered here. We assume that sequential ionization proceeds from the weakest to the strongest bound electron. The first electron (e 1 ) is ionized from the highest occupied molecular orbital (HOMO) with ionization potential I p1 followed by the second (e 2 ) from HOMO+ with ionization potential I p2 , once the first electron has been born. Simultaneous ionization is not investigated here. Atomic units are used unless otherwise noted.
The unperturbed cluster electrons are assumed to be bound by a short-range, spherical well potential. The potential of the ionizing, field dressed cluster is approximated by using the electrostatic potential of a metallic sphere in the presence of laser fieldẑF(t) and already ionized electrons [22]. The potential V 1 experienced by e 1 tunneling from the HOMO is given by The second term in V f arises from the cluster polarization which counteracts the laser field to make the potential at the surface zero. Further, a is the cluster radius, x i (t) is the position of electron i = 1, 2 with respect to the cluster centre, and r i = ∥x i ∥.
The second electron tunnels from HOMO+ to the continuum through the two-electron potential The first and second term in equation (2b) represent the combined potential of e 1 and of its image charge on e 2 . The last term ensures that the surface charge remains zero, as is required for an isolated sphere [22]. The position of e 1 in V 12 (x 1 , x 2 ) modifies the tunneling barrier for e 2 and results in CB.
Note that each electron during tunneling leaves an attractive ionic (self-image) charge behind. This results in additional pre-exponential corrections to the tunneling rates, which are not considered here. We focus on a relative comparison of tunnel ionization with and without CB, where such effects are immaterial. However, they could be included following the methods developed in atomic ionization theory [23,24].
As rigorous two electron ionization calculations of complex systems are out of reach, the wavefunction needs to be approximated. The required simplification emerges from the adiabatic nature of tunnel ionization. Tunneling theory rests on the fact that the under-the-barrier wavefunction is established practically instantaneously [28,29]. Therefore, e 2 can be assumed to respond instantaneously to e 1 by finding a new quasi-equilibrium wavefunction under the barrier, resulting in a modified tunnel ionization rate. In contrast, the effect of e 2 on the continuum wave function of e 1 is negligible compared to that of the laser. This can be understood as follows. The magnitude of the force of the two electrons on each other has to be equal; however their response can be quite different according to the different eigenfunction spectrum in which they are embedded. While e 1 experiences the classically allowed continuum and approximately follows the laws of classical physics, e 2 is composed of non-classical under-the-barrier states; the tunneling wavefunction responds practically instantaneously to changes of the tunneling barrier. As a result, x 1 in V 2 can be treated as a variable instead of a coordinate; the two electron problem is broken down into a series of one electron problems with e 1 occupying all allowed continuum positions with a probability to be determined below.
We continue with the framework for calculating tunnel ionization of e 1 and e 2 . In the quasi-static limit the potentials V i are frozen and 3D Wentzel-Kramers-Brilluoin (WKB) theory is applicable [25,26]. We use a slightly modified 3D WKB approach that builds on the slowly evolving wave approximation, a concept used in non-linear optical pulse propagation [27]. A detailed derivation and discussion of this approach can be found in the Supplementary Material 3 (stacks.iop.org/JPPHOTON/2/034007/mmedia). The tunnel ionizing wavefunction is found to be where in the following discussion the index i = 1, 2 refers to e 1 and e 2 . The pre-exponential factor in equation (3) is obtained from matching tunneling wavefunction with the initial, field free ground state. For the sake of simplicity we assume l = 0 angular momentum HOMO and HOMO+ ground states, A i exp(−κ i r i ), where A i is a normalization constant and κ i = 2I pi . The classical action S 0i and tunneling wavevector k zi (i = 1, 2) are determined by Equation (4) is solved numerically by the method of lines and Runge-Kutta integration [30]. Integration is performed subject to the initial condition S 0i (z i = a) = −κ i x 2 i + y 2 i + a 2 determined by the exponent of the field-free ground state. From matching decaying and oscillating parts of equation (3), the tunneling electron current density j zi for potential V i (i = 1, 2) is found to be where ℜ denotes the real part and x mi = (x i , y i , z mi ) represents the transverse plane at the outermost turning point, which is defined by where w i (z mi ) =˜j zi dx i dy i is the ionization rate. Expressed in words, ionization is considered to be finalized beyond some transverse plane at z mi where the wavefunction is oscillating at all points x i , y i . Finally, the dependence of j z2 on e 1 will be specified below. Next, we put together all the elements developed so far and discuss our complete theoretical approach; for a summary see schematic in figure 1. First, e 1 tunnels through the potential given by equation (1) and a transverse wavepacket ψ 1 (x m1 , t b ) is born in the continuum at time t b (equation (3)). The ionized wavepacket is driven by potential (1) and evolves into ψ 1 y 1 , z e1 (t b , t)). During continuum propagation the time dependence of V 1 needs to be accounted for. Propagation is done quasi-classically in two steps. The longitudinal evolution z e1 (t b , t) is determined from solving the classical equation of motionz e1 (t) = F 1 (z e1 , t). Here, F 1 (z 1 ) = −∂ z1 V 1 ((x 1 , y 1 ) = 0, z 1 ) and the double dot denotes the second time derivative. The transverse component of the second term in equation (1) yields higher order contributions and is therefore neglected, in addition to the pre-exponential transverse dependence of the initial wavepacket at z m1 . Taylor expansion of the classical action about the z-axis yields, to second order, Transverse quantum spreading of e 1 proceeds as for a free Gaussian wavepacket. Applying the analytical expression for quantum diffusion of a Gaussian wavepacket [31] to the approximated ψ 1 (x m1 ) yields the transverse wavefunction ψ 1 (x e1 (t b , t)), from which the current density is found to be From equation (7) follows the rate at which e 1 , after having been ionized at time t b , passes through the transverse plane at z e1 (t b , t) at time t: The probability for e 2 to experience tunneling potential The solution of equation (4) for this potential, together with equation (5) for i = 2, determines tunneling current  e1 around xe1(t b , t) is jz1 (xe1(t b , t))dx1dy1dt b ; image charges (red dots) + charge of e1 (blue dot) modify the potential e2 sees from (t b , t), x2) and move the turning point of e2 from z (n) . The resulting conditional current density jz2 (xm2|xe1(t b , t)) of e2 is reduced, as shown here for e1 at xe1(t b , t) = (0, 0, ze1(t b , t)). density j z2 (x m2 |x e1 (t b , t)), from which the ionization rate is obtained via w 2 (x m2 |x e1 (t b , t)) =˜j z2 dx 2 dy 2 . Here, the vertical bar denotes conditional probability; w 2 is the tunneling rate at which e 2 is emitted into the transverse plane at z m2 at time t, conditional to e 1 being at x e1 (t b , t). Combining the ionization rates for e 1 and e 2 yields the sequential two electron ionization yield as Here, t l and t u are lower and upper time limits; in the case of a laser pulse they are ±∞. However, as CB is important only when both electrons are born within one laser half-cycle, it is also instructive to investigate half-cycle intervals, where t l , t u represent the times of two adjacent field nodal points. In the absence of CB, V 2 = V f , ionization of e 2 becomes independent of e 1 and the sequential double ionization yield transitions into Note that there is a probability for e 2 to ionize from HOMO+ before e 1 ionizes from HOMO. While in most systems it is small, in some systems it can be significant. This channel can be accounted for by flipping the roles of e 1 and e 2 in equations (8) and (9), but we do not consider it here for the sake of simplicity. The remainder of this discussion focuses on the numerical results. Ionization potentials for small clusters were estimated using a spherical drop model applied to tantalum clusters, where I p = Φ + A/r, Φ is the work function, r is the cluster radius, and A is an empirically determined parameter [33]. As CB manifests predominantly in the exponent of the tunneling current, pre-exponential contributions are not considered in the evaluation of the above formalism.
In figure 2, a cluster with radius a = 6 is exposed to a field F 0 sin(ωt) over a half-cycle 0 ≤ t ≤ T 0 /2 with ω 0 = 0.057, T 0 = 2π/ω 0 , and F 0 = 0.03. Additionally, birth time and coordinates of e 1 are t b = 0.2T 0 and (x 1 , y 1 ) = 0. In figure 2(a) the longitudinal evolution of e 1 is plotted versus time t. In figures 2(b-e) ℜ{S 02 (x 2 , 0, z m2 )} is plotted versus x 2 for increasing times t marked by the empty (red) circles in figure 2(a); the classical action is connected to ionization rate via equation (5). The full (blue) and dashed (red) lines represent the results in the absence and presence of CB, respectively. The magnitude of ℜS 02 drops with increasing t, as the electric field decreases. CB bores a hole in the centre of the classical action, where ionization is otherwise maximum. As a result, ionization suppression is strongest when e 1 is born on-centre. For e 1 born off-centre, CB is less effective, as the hole is shifted off-centre where ionization is already less likely. Further, the radial symmetry of S 02 is broken. For the parameters chosen, CB remains pronounced throughout the whole time. This appears at first counter intuitive, as e 1 is moving away from the cluster with increasing t. However, the resulting drop is compensated by the fact that the tunneling barrier of e 2 is also shifted outwards due to a drop in field strength.
In figure 3(a) the importance of CB as a function of wavelength and cluster radius is investigated for F 0 = 0.02. This is done by calculating sequential two electron ionization yields n cb 12 and n n 12 defined in equations (8) and (9) over a half cycle. The ratio ∆n 12 = n n 12 /n cb 12 is plotted versus λ, where curves 1-4 refer respectively to cluster radii a = 6, 12, 25, 100. Measuring the ionization yield over one half cycle gives an  upper limit for CB, which can be approached with single cycle pulses, where ionization occurs predominantly during one half cycle. Figure 3(a) shows that CB is most pronounced in the short wavelength and small cluster limit. This is intuitive, as electrons have less opportunity to avoid each other in smaller geometries. Moreover, the velocity of e 1 scales as F 0 /ω 0 , causing e 1 to be driven more quickly from the parent cluster at longer wavelengths, giving it less time to block e 2 . Note however, that the cluster size trend does not extend down to atomic sized systems. In atoms the difference between I p1 and I p2 is much larger, making it unlikely that both electrons can tunnel during the same half cycle. Conversely, as cluster size increases, ionization potentials decrease and so too does the relative difference between them, yielding a progressively weaker CB effect due to similar ionization potentials and the aforementioned capacity for spatial avoidance. Finally, for λ 0 = 800nm and for a = 6 the influence of CB is substantial, reducing n 12 by more than two orders of magnitude.
In figure 3(b) ∆n 12 over one half cycle is plotted as a function of F 0 for ω 0 = 0.057 and a = 6. The importance of CB increases with decreasing field strength, for which e 1 stays longer in the vicinity of the cluster, thus blocking ionization of e 2 more effectively.
Finally, in figure 4 total ionization n 12 over a Gaussian laser pulse, F(t) = F 0 exp(−(t/τ ) 2 ) cos(ω 0 t), is plotted as function of pulse duration τ , with laser parameters F 0 = 0.03 and ω 0 = 0.057, and cluster radius a = 6. Full (blue) and dashed (red) lines refer to ionization in the absence and presence of CB. The green empty circles indicate a τ 2 -fit which is expected for sequential ionization n n 12 . We find that for τ /T 0 < 2 the effect of CB becomes very pronounced. It results in a sharp drop in ionization below the τ 2 -fit, whereas in the absence of CB the drop in ionization begins to level out. For single cycle pulses ionization occurs dominantly during one half cycle, where the effect of CB is most pronounced. In longer pulses, e 1 and e 2 can be born in different half cycles, in which case the effect of CB is negligible. As a result, the effect of CB measured over a multi-cycle pulse becomes diluted.
We have shown that Coulomb blocking suppresses sequential two electron tunnel ionization; in sub-nm clusters driven by single-cycle, near-IR laser fields, suppression by 2-3 orders of magnitude was found. Coulomb blocking is not limited to two electrons. If more than two electrons can tunnel during a half-cycle, each electron will suppress tunneling of each subsequent electron. Indeed, our model is in principle extensible to include more than two active electrons. It should be noted, however, that each additional electron significantly increases computation time. Our analysis revealed various ways to experimentally identify Coulomb blocking. First, sequential two electron tunnel ionization was found to be suppressed strongly for pulse durations decreasing towards the single cycle pulse limit. Second, tunneling of e 2 is strongly suppressed along the direction in which e 1 has tunneled, which could be measured in COLTRIMS (cold target recoil ion momentum spectroscopy) experiments [32]. Finally, the amount of suppression depends on the birth times of electrons 1 and 2 varying over a laser half-cycle resulting in a modulation of two-electron emission on an attosecond time scale. Such signatures could be resolved by two-colour attosecond experiments.
The authors wish to acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). I. Schubert acknowledges the German Academic Scholarship Foundation for financial support. Some simulations were performed on the Mammouth Parallèle cluster provided by Compute Canada and maintained by Calcul Québec at l'Université de Sherbrooke.