Quantum battles in attoscience: tunnelling

Abstract What is the nature of tunnelling? This yet unanswered question is as pertinent today as it was at the dawn of quantum mechanics. This article presents a cross section of current perspectives on the interpretation, computational modelling, and numerical investigation of tunnelling processes in attosecond physics as debated in the Quantum Battles in Attoscience virtual workshop 2020. Graphic abstract


Introduction
The discovery of the quantum tunnelling phenomenon almost 100 years ago has not only opened up many new avenues and applications. It has also kept quantum physics researchers busy since then, trying to define the temporal resolution of the process [1,2]. Early experiments were focused on photons tunnelling through potential barriers, such as Ref. [3] for example. But with the advent of attosecond science [4] the question "Does tunnelling take time, and if yes, how much? " has gained a lot of new interest, since electron dynamics often include quantum tunnelling portions, be that in biological processes such as photosynthesis [5] or charge transport in semiconductors [6], tunnelling ionisation as the first step for high-order harmonic generation (HHG) spectroscopy [7], photoelectron holography [8], laser induced electron diffraction (LIED) [9] or many more.
The temporal resolution of quantum tunnelling is still heavily debated [10][11][12][13] and thus presented an interesting topic for a debate at the Quantum Battles in Attoscience virtual workshop 2020 [14]. The aim of the Battle sessions was "an open debate on a contentious topic involving several early career researchers ('combatants') and the entire audience of attendees" [15]. To that effect, the combatants prepared a scaffolding structure of the debate on "tunnelling", defining three main topics: (a) Physical observables and typical experiments (presented in Sect. 2 of this article), (b) Nature of Tunnelling (see Sect. 3), and (c) Theoretical approaches to quantum tunnelling time (in Sect. 4). Each topic was introduced with an overview presentation, followed by a a e-mail: c.hofmann@ucl.ac.uk (corresponding author) free debate among all combatants, moderated by Prof. Jonathan Tennyson, UCL, and included both questions among the combatants as well as live audience questions. The result was a highly interactive and lively debate [16].
This perspective article offers a text-form of the live debate [17], supplemented with additional references and explanations.

Physical observables and typical experiments
The guiding questions for this first topic are: What are physical observables, typical measurements, and what are the characteristic physical systems under investigation? What other aspects of these particular systems influence the interpretation of tunnelling time studies?

Overview
When it comes to experiments investigating the temporal resolution of a quantum tunnelling particle, there are typically two kinds of experiments: (a) Bose-Einstein-Condensates (BEC) of atoms trapped in optical lattices, with various manipulations on them to measure tunnelling from one lattice site to the next [13,18]. Since the particles in question are entire atoms, their temporal resolution for the dynamics is in the range of microseconds. And (b) attosecond angular streaking (also known as attoclock) type [19][20][21] exper- iments, a technique developed in strong-field attosecond physics, where electrons tunnel ionise from a bound state through the potential barrier which is created by the interaction of the strong laser field with the binding Coulomb potential of atoms. These are on the attosecond regime since electrons are tunnelling, and the main focus of the here following debate.
On a fundamental level, what we are interested in is the temporal resolution of a wave packet hitting a potential barrier, and then a part of that wave packet tunnelling through, such as schematically illustrated in Fig. 1. However, this exactly creates several challenges in trying to time this process compared to other timings of wave packets, such as for example group delay in photonics. The peak of the wave packet is not conserved, since the incoming (or bound state) wave packet is split into a reflected and a transmitted part. The potential barrier essentially acts as an energy-dependent filter, such that the spectra of the two resulting wave packets are significantly different [10]. A wave packet also always corresponds to a probability distribution, and in consequence it is difficult to define a clear starting and ending (or entrance and exit) point, more on that in Sect. 3. Furthermore, in strong-field attoscience scenarios we are tunnel ionising from a bound state, where of course parts of the wave function even in its field-free ground state always are "under" the barrier, without any tunnelling occurring. Additionally, approaches such as Wigner-like, scattering and resonance phase times [22] which are commonly applied to single-photon ionisation [23] are not applicable either, again because of the chirped propagation of the electron wave packet and the energy filtering of the potential barrier. While we are on the topic of potential barriers, it is also worthwhile noting that this classical picture of the potential barrier only emerges if the laser field is treated in the length gauge [24][25][26][27].
The physical observable for measurements (and calculations often, too) of strong-field tunnel ionisation are momentum distributions of photoelectrons [20,21,28] or momenta of atoms [13,18]. Momentum is of course a standard quantum mechanical observable corresponding to a unitary operator, whereas time itself is a parameter of the Schrödinger equation and thus not an observable as such. Therefore, a relation between measured (or calculated) momenta and the timing of the tunnelling process needs to be established through theoretical understanding of the quantum tunnelling process.
In the experiment by Fortun and co-workers a rubidium BEC is oscillating in an optical lattice. In a pumpprobe-type approach, the lattice is turned off at different intervals after the initiation of the oscillation and the instantaneous momentum of the atoms carried them flying towards a position-sensitive detector. The tunnelled wave packets appeared delayed with respect to the reflected wave packets [18]. In the experiment by Ramos and co-workers, a quantum simulation of the Larmor clock [29][30][31], one of the well-known theoretical approaches to predicting the tunnelling time [32], was realised causing precession of the spin of the rubidium atoms while traversing a potential barrier. This spin precession was then mapped onto different states according to the angle of rotation and separated by a Stern-Gerlach measurement [13].
In attoscience experiments utilising the attoclock method [33,34], the rotation of the nearly circularly polarised vector potential A mimics the hand of a clock. The path of a photoelectron after tunnel ionisation is dominated by the interaction with the laser field [35], and thus neglecting all other corrections and perturbations, the final asymptotic momentum p f is determined by the vector potential at the time when it first exits the potential barrier and enters the continuum t 0 , through the conservation of canonical momentum where p 0 denotes a possible initial momentum. Hence, the final momentum angle acts as a clock for the exit moment in time. However, this angle to time mapping is subject to several corrections, some easy to describe and include in calculations, others more elusive to quantify and thus the topic of ongoing research. A non-exhaustive list of corrections, approximations, and other issues include: the Coulomb force of the parent ion induces an angular shift [11,21,36,37]; the ellipticity, pulse envelope, pulse duration, and carrier-envelopeoffset phase are wave form parameters which affect the photoelectron trajectories; and depletion mixes in with pulse duration and the intensity of the applied field [38,39] for topics which mostly have been dealt with in great detail and comparable results; the experiment does not have access to any "start" signal of the tunnelling process, only the exit point [10,28]; nonadiabatic effects influence the ionisation rate, energy at tunnel exit, initial momentum p 0 distribution, and the location of said tunnel exit itself [39][40][41][42][43]; multi-electroneffects are ignored in most calculations [44][45][46][47][48]; models including non-classical characteristics of the trajectory which can be compared against experimental data are still being developed [20,38,49]; and the orbital angular momentum of the bound state has an effect on the strong-field ionisation [50][51][52][53] for issues which are more elusive (although this categorisation is not definite). There is still a lot of work necessary to properly disentangle the different contributions which lead to various angular shifts, sketched in Fig. 2 of the measured Photoelectron Momentum Distribution (PMD), until we can Illustration of photoelectron momentum distribution for ellipticity 0.87, clockwise helicity, projected to the plane of polarisation. Single Classical Trajectory (SCT) models assuming instantaneous tunnelling predict an angle offset away from the pure −A(tmax), but the measured angle offset might be even larger than that. Adapted from [10] Fig. 3 Time-dependent Schrödinger equation (TDSE) calculation of photoelectron momentum distributions for hydrogen ionisation. Left: idealised attoclock with a single cycle pulse and circular polarisation leads to a unique final momentum probability distribution peak. Pulse duration ≈ 1.6 fs FWHM, peak intensity 0.86 × 10 14 W/cm 2 , wavelength 800 nm, with clockwise helicity. Right: A multicycle pulse yields two main blobs with Above Threshold Ionisation (ATI) rings from the inter-cycle-interference. Pulse duration ≈ 6 fs FWHM, peak intensity 1.5 × 10 14 W/cm 2 , wavelength 770 nm, with clockwise helicity. Adapted from [11,54] be sure of the remaining angle offset and it's relation to tunnelling time.

Debate
- Figure 3 exemplifies the pulse duration and wave form dependence, as well as an energy dependence between the different ATI rings in the long pulse case, which show different angular maxima [55][56][57]. This raises concerns about the validity of one single time (rather than a distribution of times) extracted typically from data, thus averaging over the energy dependence. In attoclock experiments, the carrier-envelope-offset phase (CEP) was not stabilised [20,21,28]. Additionally, the orientation of the polarisation ellipse in the lab frame was chosen such that the observable of interest (angular shift mostly parallel to the major axis of polarisation) was orthogonal to the direction with the biggest experimental noise (along the gas jet direction, thus chosen for the minor axis of polarisation) [58]. Both of these effects wash out ATI interference. For the CEP influence in particular, the interplay between ellipticity and pulse duration is critical. For the largest field strength to be following the polarisation ellipse (desired in attoclock experiments [10]) rather than the CEP [33,34], the pulse envelope must be long enough relative to the ellipticity reducing the field strength within a quarter cycle. Furthermore, ATI rings result from interference created by many laser optical cycles, highlighting the difficulty of defining a "single time", or even relative time intervals with respect to local maxima of the field strength or other possible references. -Most often, tunnelling time calculations tend to use only a single peak point in the momentum distribution [12,38,59]. However, based on this discussion it would seem more appropriate to extract the tunnelling time from the full momentum distribution, which contains much more information regarding the tunnelling process [32,39]. We note that such work has been carried out in a recent publication [60], which assesses the whole momentum distribution instead of just a single offset angle. -An audience question is brought in: How is the peak of the PMD determined precisely, since the maximum in a 2D distribution is not the same as the maximum in the 1D angular distribution? Of course the strictly linear angle-time relationship is only exact for circular polarisation. For any other polarisation, the elliptical geometry introduces corrections and needs to be taken into consideration [61]. These effects as well as the influence of integrating over the radial component in the 2D distribution were double-checked against. The resulting shifts in the extracted values were smaller than or of the order of the reported error bars for experimental data. Nevertheless, it is important to keep in mind that different coordinate system transforms and peak angle extraction methods lead to significant shifts in the extracted angle and thus the interpreted delay time [62]. Regarding the third component, the laser propagation direction, so far no significant difference has been found between a projection to or a cut along the polarisation plane. Of course this requires that no extra physics becomes important along this third component, for example the influence of the magnetic field must be negligible [26,63]. Of course, an energy-resolved angular distribution would avoid the integration over at least one of the components and thus make the peak search less dependent on geometry and coordinate choices.

Nature of tunnelling
The guiding questions for the second topic are: What is the nature of tunnelling at the classical/quantum intersection? What is the "beginning" and "end" of tunnelling, and how do we define it? What are classical or quantum trajectories?

Overview
Quantum tunnelling is a wave phenomenon, and the time-dependent Schrödinger equation (TDSE) is an equation for the probability amplitude wave (wave function). But this description makes it difficult to define where and when tunnelling exactly starts. Tunnelling itself is natural in quantum mechanics, it is only when we look at it from a classical perspective that there is a "forbidden" region in the potential barrier. In the classical domain, trajectories are well defined in space and time, but can they tunnel? The semiclassical models typically use classical trajectories to describe the motion of an electron after it has been released from an atom, usually by tunnelling ionisation. Is a synthesis of these two worlds like the sketch in Fig. 4, aiming to retain the quantum physics behaviour with the clarity of trajectories, possible? It is clear that such a synthesis is not a simple task. Indeed, in order to calculate the classical trajectory, i.e., to integrate Newton's equation of motion, both starting point and the initial velocity are needed. However, Heisenberg's uncertainty principle imposes a fundamental limit to the accuracy with which the values of the position and momentum, as well as of any other canonically conjugate variables, can be simultaneously determined. Nev-ertheless, the application of the quasiprobability distribution allows to obtain information about both the position and momentum from the wave function. The most widely known examples of quasiprobability distributions are the Wigner function and Husimi distribution. We note that the Wigner function has already been used for description of strong-field processes, see, e.g., Refs. [64][65][66]. However, to the best of our knowledge, the Wigner function has not yet been applied to the combination of the quantum and trajectory-based description in strong-field tunnel ionisation, although a similar method has been proposed for the case of attosecond pulse single-photon ionisation with subsequent streaking of the photoelectron wave function [67]. A recent and successful attempt of such combination was made in Ref. [49] using Gabor transform.
Doing so still begs the question, which quantity best characterises the onset of tunnelling?

Debate
-Quantum particle description is necessary for tunnelling to occur in the first place, but the potential barrier defines local properties which are significant for classical systems. So we need a combination of both, quantum tunnelling feature with the classical flavour of understanding if we aim for any kind of temporal resolution of a tunnelling process. The challenge is then to find one single picture for the entire process. Instead of relying on real-space trajectories, including complex space and time enables the tunnelling phenomenon, resulting in a quantum trajectory with clearly defined entry and exit to the barrier [68], as well as corresponding times (more on this method in Sect. 4).
On the other hand, measurements can always only find real observables, thus fully complex calculations must find their way to the real axis somehow, where the propagation of a photoelectron wave packet is very well described by classical methods [69]. But since the experimental observables typically are momenta, purely quantum models which operate in complex space and time can still be used for the purpose of comparison, as long as they can predict a final momentum distribution. -One huge assumption in experimental approaches based on the "attoclock" principle is the "starting time", relative to which the tunnelling delay is calculated. This is typically chosen to be the maximum of the electric field, since that moment corresponds to the highest probability of tunnelling [10]. However, this assumption might be missing out on half of the effect [41], and the tunnelling process in strong-field ionisation might be a symmetric problem relative to the (local) field maximum [70]. Publications which attempted to identify a physical starting point have found other values, typically before the maximum is reached [41,42,71,72]. -If we consider fully quantum models which describe both the tunnelling transition from bound to ionised state and the propagation afterwards in one, it becomes important to distinguish tunnelling from over-the-barrier (OBI) ionisation. In experimental approaches it is generally not possible to aposteriori separate these two contributions to the total momentum distribution, and the same limitation is also true for numerical solutions of the TDSE [73]. However, in theoretical calculations based on trajectories, these two processes can easily be differentiated. The semiclassical models naturally distinguish between the tunnelling through a potential barrier and the over-barrier-ionisation. Indeed, when the field strength is so high that the potential barrier formed by the laser field and the ionic potential is suppressed, it is impossible to find the starting point of the electron trajectory using field direction model (see, e.g., Refs. [74][75][76][77]) or the separation of the static tunnelling problem in parabolic coordinates [78]. In this case it is usually assumed that the electron starts at the top of the suppressed potential barrier, and the difference between the ionisation potential and the energy at the top of the barrier ΔE = −I p − V max is transferred to the initial longitudinal velocity of the departing electron: Non-adiabatic effects, i.e. effects beyond the quasistatic approximation which are due to the timedependent changes in the strong field, also play into these definitions. For example, at which energy or distance can a photoelectron exit the potential barrier [43,79] or when does the onset of OBI occur?
The so-called backpropagation method [39,42,43] is one hybrid approach which utilises the full quantum power of the TDSE for the tunnel ionisation but then retroactively adds the power of classical trajectories to also distinguish between OBI and tunnelling (more on this method in Sect. 4). -An audience member suggests that localised position measurements would be able to distinguish tunnelling and OBI, since only OBI would be detectable. This gedanken experiment however would require a detector positioned at the atomic potential barrier, which is unfeasible in any kind of experimental setup since detectors require some time-of-flight information and are placed a significant distance away from the interaction region, of the order of several centimetres at least [58,80]. There have been some theoretical studies using virtual detectors in combination with TDSE solutions [71,72], but those again can not distinguish of course. This is because both under-the-barrier and over-the-barrier transmission causes a probability flux of the wave function, which a hypothetical detector would be able to pick up without being able to distinguish between those two types of transmission. This remains the case also in different tunnelling scenarios such as in a tunnelling junction where a macroscopic, position resolving detector might be feasible. -The last point for this topic is concerning representation of a quantum wave function by using trajectories. Fundamentally, we are trying to study the behaviour of a wave function doing something interesting. In the most simple trajectory approach, the entire wave function is represented by a single simplified wave packet (i.e. a Gaussian) with the associated trajectory of its wave packet peak mapping the motion of the expectation value of the wave function (similar to a group velocity approach) [81]. However, this approach can not describe a wave packet being split into a reflected and a transmitted part, and thus would either always remain bound or the bound state is fully depleted. Trajectories representing a skewed wave packet [82] would present a more generalised version. Even more accurate are descriptions that employ a large ensemble of trajectories following the probability distribution of the underlying wave function [83]. The question is then: Will an ensemble of (classical or quantum) trajectories not only represent an instantaneous probability distribution derived from a wave function, but also its dynamics over time?
This question was addressed in Ref. [84] that studies the validity of the two-step semiclassical model disregarding quantum interference but accounting for the Coulomb field for strong-field ionisation. The Ehrenfest theorem [69] (see, e.g., Ref. [85] for a textbook treatment), which establishes quantum mechanical analogues of classical Hamiltonian equations, was applied in Ref. [84]. Furthermore, the analysis of Ref. [84] is based on a quantitative comparison of the electron momentum distributions obtained within the two-step model and by numerical solution of the TDSE. Reference [84] introduces the measure for the deviation of the dynamics of an ensemble of classical trajectories from the Ehrenfest's theorem. This measure is the relative deviation between the force at the average position of the ensemble of trajectories and the average of the forces on the ensemble. A correlation was found between the invalidity of the two-step model and the deviation of the dynamics from the Ehrenfest's theorem. The general trends for the applicability of the twostep model in terms of laser intensity, wavelength, ellipticity, as well as in terms of the potential properties are identified in [84]. However, this study is done in the two-dimensional (2D) case and needs to be extended to the 3D one.
What are theoretical approaches used to investigate quantum tunnelling times? What are their various characteristics, advantages and disadvantages?

Overview
For the overview, a brief and non-exhaustive list of different calculation approaches to the tunnelling time are given. They are categorised with regards to their theoretical foundation.

Quantum methods based on time-dependent Schrödinger equation
First are numerical solutions to TDSE. A common advantage of all these methods is that they are fully quantum calculations for the entire process. Further individual characteristics, advantages and disadvantages can be summarised as follows. TDSE calculations which employ Coulomb vs Yukawa potentials [21,38] found that attoclock signal shows a prominent offset angle with Coulomb binding potential, while the offset angle vanishes for a Yukawa potential. This comparison offered an indirect proof of instantaneous tunnelling by comparing the results depending on the two different binding potential of the parent ion.
The numerical saddle-point method [62] uses a trajectory-free language and establishes a connection between the final momentum of the photoelectron and the numerical saddle-point time for the full Hamiltonian including the Coulomb potential. It supports the conclusion of instantaneous tunnelling. However, this method is gauge dependent.
The functional derivative method [70] investigates the instantaneous ionisation probability as a functional derivative of the total ionisation with respect to the wave form of the ionising field, but does not map directly to any experimental observables. It is gauge independent, and found vanishing delay (or vanishing delay asymmetry with respect to the local peak in the field).
Bohmian mechanics [86] present a mapping from the quantum world to the trajectory language. However, the calculation is guided by a pilot wave not pertaining solely to the (eventually) ionised part of the wave packet near the tunnel exit, thus potentially giving false tunnelling information. A separation of the (eventually) ionised part and bound part of the wave packet near the tunnel exit is, unfortunately, impossible, due to quantum nonlocality.

Quantum methods based on strong-field approximation
Strong-field-approximation (SFA) [87][88][89] based quantum methods describe ionisation as a transition from an initial state unaffected by the laser field to a Volkov state, i.e., the free electron wave function in an electromagnetic field. Therefore, the SFA disregards the  [90] intermediate bound states and the ionic potential (e.g., Coulomb interaction) in the final state. Presently, several SFA-based quantum approaches are developed. Typically, these approaches decide which force dominates the trajectory of a photoelectron based on its position in space and use the corresponding approximations. This separation and reduction of the acting forces allows for analytic calculations. The imaginary part of the saddle-point time in SFA calculations relates to the inverse tunnelling rate, while the real part in these models is often taken as the tunnel exit time.
The analytic R-Matrix (ARM) method [36,38] separates space into an inner region (Coulomb & Laser field considered) and an outer region (Coulomb field neglected, eikonal-Volkov approximation), as illustrated in Fig. 5. The disadvantage of this method is the challenge of choosing proper integration contours for each trajectory.
The under-barrier recollision theory [90] specifically includes interference between under-barrier rescattered and direct trajectories, as shown in Fig. 6. This leads to a shift in the momentum wave packet peak, which can be interpreted as a delay. However, this method ignores Coulomb corrections.

Hybrid quantum-classical method
The backpropagation method [39,42,43] is a hybrid quantum-classical approach offering a unique perspective on the tunnelling process. It combines a fully quantum calculation of the ionisation process with for- Fig. 7 Concept of the backpropagation method ward propagation utilising TDSE solution, followed by a transcription of the resulting ionised quantum wave packet into classical trajectories, and a subsequent propagation of the trajectories backward in time, see Fig. 7 for a sketch. Another variant of the backpropagation method would be putting a sphere of virtual detectors [91][92][93][94][95][96] around the target, where the flux is converted into classical trajectories during the laser pulse on the fly [53,97].
Why backpropagation? Firstly, as everyone agrees, tunnelling is a purely quantum process. Introducing a tunnelling barrier into the description of tunnelling ionisation, however, brings in clearly classical elements into the picture. Namely, tunnelling is now depicted with local tunnelling exit positions and momenta, which calls for a classical formulation. Secondly, due to quantum nonlocality, the portion of the wave packet that would eventually be freed and the portion that would finally remain bound can not be separated during the tunnelling process. A separation is only possible in the far field, when these two portions are spatially detached. These are exactly the design philosophy of the backpropagation method, a hybridisation of quantum forward and classical backward propagation. It combines the advantages of the quantum and classical methods by offering the capability to include the full Hamiltonian and quantum tunnelling dynamics while retaining the local information from the classical trajectories. It also naturally includes nonadiabatic tunnelling effects, automatically remove the offset angle from Coulomb effects, and retrieves the electron characteristics at the tunnel exit.
The classical backpropagating trajectories may be stopped whenever a certain condition is met, which defines the tunnel exit, yielding highly differential information of the tunnel exit. In this manner, the backpropagation method may act as a common ground to compare different definitions of tunnelling. It was found that a vanishing tunnelling time results if the tunnel exit is defined in the momentum space when the velocity of the trajectory vanishes in the instantaneous field direction (the velocity criterion), while defining the tunnel exit as a certain position in the coordinate space (the position criterion) gives rise to a finite tunnelling time [39,43]. Different definitions of the tunnel exit were thus believed to be the origin of the tunnelling time debate. It was further argued that the position criterion leads to inconsistencies and difficulties and thus the velocity criterion is favoured as the definition of the tunnel exit, and the tunnelling time delay should thus vanish [39,43].
The backpropagation method has further enabled a study of the tunnelling time delay induced by orbital deformation [53] and a subcycle time resolution of the linear laser momentum transfer, where a coupling between the nondipole and nonadibatic tunnelling effects was found [97].

Semiclassical methods
Semiclassical methods apply classical trajectories to describe the motion of an electron after it has been released from an atom or molecule by the laser pulse. The two-step [98][99][100] and the three-step [101,102] models are the most widely known examples of the semiclassical approaches. These models do not account for the effect of the ionic potential on the electron motion in the continuum. Presently there are many trajectory-based models that do account for the ionic potential in the classical equations of motion. Among these are: Trajectory-based Coulomb SFA (TCSFA) [103,104], Quantum trajectory Monte-Carlo method (QTMC) [105], Coulomb quantum orbit strong-field approximation (CQSFA) [106][107][108][109][110][111], semiclassical twostep model (SCTS) [112], Quasistatic Wigner method [20], etc. The three-step model using complex classical trajectories [68] and the classical Keldysh-Rutherford model [37] are closely related to this group of models.
Using a purely classical description of the electron motion it is not possible to describe the quantum interference effect in the photoelectron momentum distributions and energy spectra. Recently substantial progress has been achieved along these lines. Along with some other approaches, the TCSFA, QTMC, CQSFA, and SCTS models account for interference effects. In these approaches every classical trajectory is assigned to a certain phase, and the contributions of different trajectories leading to a given final electron momentum are added coherently.
The TCSFA extends the well-known Coulombcorrected strong-field approximation (CCSFA) [113,114] by treating the laser field and the Coulomb force acting on the electron from the ion on an equal footing. The TCSFA accounts for the Coulomb potential in the phase of every trajectory within the semiclassical perturbation theory. The same approach is used in the QTMC model. In contrast to this, the SCTS and the CQSFA models account for the Coulomb potential beyond the semiclassical perturbation theory.
The quasistatic Wigner method [20] employs the concept of the dominant quantum path. Using the spacetime propagator, the quasistatic Wigner method considers the propagation of the electron wave function that originates from the initial bound state in the classically forbidden domain. The quasistatic description of the laser field is used in Ref. [20]. The phase of the quantum mechanical propagator determines the most dominant path along the tunnel channel, and therefore, determines the Wigner trajectory. The Wigner trajectory is merged with the corresponding classical trajectory in the continuum, see Ref. [20]. In this way the  [68] quasistatic Wigner method determines the initial conditions for the classical trajectory. It should be emphasised that the initial conditions include not only an initial momentum, but also a time delay. However, this method reduces the wave packet to a single trajectory. It should also be noted that the Wigner time is illdefined in the tunnelling process [10].
Since real-valued trajectories are not able to describe tunnel ionisation, the complex-time-and-space model [68] employs complex trajectories, as illustrated in Fig. 8. This approach was applied to the HHG process in Ref. [68]. All components of the three-step model are described in Ref. [68] within a single consistent trajectory framework. The trajectories are sampled from an initial Coulomb eigenstate, and the time propagation is performed using the final value coherent state propagator (see Ref. [115]). As a result, the model provides a unified and seamless trajectory description of the ground state, tunnelling, and collision process. The model shows quantitative agreement with fully quantum results. However, the contour in the plane of complex time, which is necessary to implement the model, has to be chosen manually.

Classical methods
And finally, purely classical models are still also developed and used often. The Keldysh-Rutherford model [37] applies the famous Rutherford scattering formula taking the vector potential of the laser pulse as the asymptotic electron velocity and the Keldysh tunnelling width as the impact parameter. The model was tested by comparison of its predictions with the numerical solution of the TDSE using the hydrogenic potential and the screened (Yukawa) potential. In the latter case the action of the Coulomb field was gradually switched off. The striking similarity between the attoclock offset angle and the Rutherford scattering angle was revealed in Ref. [37]. The Keldysh-Rutherford model suggests that the offset angle has a largely Coulombic origin [37]. Therefore, the model is questioning the interpretation of this angle in terms of a finite tunnelling time. However, the Keldysh-Rutherford model completely neglects nonadiabatic effects, and is also limited in its validity to short pulse durations and (relatively) weak intensities which are outside the typical parameter range of experiments to date. Therefore, some further work along this direction is needed.
Classical trajectory Monte Carlo (CTMC) methods [28,45] are the classical cousin of QTMC, and often employed where interference effects are not of any key interest. Since the calculations are computationally cheap compared to TDSE solutions and fewer trajectories are needed than in QTMC to reach similar statistical quality, these methods are able to fully include the ion Coulomb potential together with the laser field during the propagation after the tunnel exit, as well as various non-adiabatic effects [116], certain multi-electron effects [47] including Stark shift and an induced dipole in the parent ion [45].

Debate
-Regarding the complex-time-and-space method [68], the question is raised how exactly the integration contour is chosen manually. Every time a trajectory orbits the singularity at the nucleus, there is a possibility for the quantum trajectory to be emitted from the bound wave packet and leave as part of the ionisation wave packet. The exact choice for after how many orbits an ionisation event happens is done by comparing to the full quantum result, for sections of initial coordinate space. Observables are then computed from the resulting trajectories. The number of loops is a discrete choice and not a fully tunable parameter. While there is a choice to match the quantum result, each discrete choice yields significantly different results, so the agreement with TDSE calculations is not entirely by construction. The interpretation of this choice is not clear yet from a physical point of view.
When it comes to separating the eventually bound from the eventually ionised part of the wave packet, trajectories far from the core in the long time limit are considered ionised. Searching for conditions (zero momentum for example) along those eventually ionised trajectories yields two complex times, labelled as tunnel entry and exit, where the difference can be interpreted as tunnelling time, see Fig. 9a. Alternatively, the difference to the field maximum can be computed, as shown in Fig. 9b. This may be required to compare the results to experiments where tunnel entry times may not be accessible but it does ignore a significant contribution to the total tunnelling time. As Fig. 9 shows, the time required for tunnelling is nonzero in both real and imaginary components. Furthermore, a single averaged result may be insufficient to characterise the tunnelling process as the distribution of times is wide and asymmetric. Two distinct classical processes are found, and the states on an attoclock measurement is. This was dealt with in the [21] study since molecular hydrogen had to be split into atomic hydrogen. In their extended data figures & tables, it is shown how initial bound states 1s or 2s result in completely different final momentum distributions. Photoelectrons ionised from 2s have much smaller absolute momenta, their distribution shows a different structure, and the event is less likely to happen. Therefore, contributions from different initial states can be separated. -Already in Sect. 3, combined approaches which offer quantum behaviour with trajectory insight have been identified as beneficial for many strong-field (tunnelling) phenomena models.
Typically, the SCTS model requires large ensembles of classical trajectories to resolve fine interference details. These trajectories are propagated, and their final momenta are binned in cells in momentum space. This is often referred to as "shooting method' [103], although this approach has nothing to do with the shooting method for solving a boundary value problem. In contrast to the TCSFA, QTMC, and SCTS, the CQSFA method finds all the trajectories corresponding to the given final momentum. This approach is often called the solution of the "inverse problem" and it allows to bypass the necessity of large ensembles of trajectories. However, the solution of the inverse problem is a non-trivial task, and, furthermore, is generally less versatile than the "shooting method". Any trajectory-based model requires specification of initial conditions, i.e., the initial electron velocity and the starting point of the classical trajectory. Indeed, these initial conditions are needed to integrate the Newton's equations of motion. The starting point, i.e., the tunnel exit, is found using the separation of the tunnelling problem in parabolic coordinates [78]. The Stark shift of the energy level that has an effect on both the tunnel exit and ionisation probability was also taken into account in the SCTS. It is generally considered in the semiclassical models that the electron departs with zero initial velocity along the laser polarisation direction v 0, = 0 and an arbitrary initial velocity v 0,⊥ in the perpendicular direction. The ionisation times and the initial transverse velocities are distributed in accord with the static ionisation rate: with κ = 2I p . The quasistatic approximation is used in Eq. (3), i.e., the static field strength F is replaced by the instantaneous value F (t 0 ). The quasistatic approximation is used in both QTMC and the SCTS. We note that many trajectory-based models use the SFA formulas instead of Eq. (3) to distribute the initial conditions of classical trajectories, see, e.g., Refs.
[103, [116][117][118][119]. This allows to investigate nonadiabatic effects in above-threshold ionisation and often leads to a better agreement with the numerical solution of the TDSE. Recently, the SFA-based formulas as distributions of the initial conditions have been validated in a systematic way [60]. It is found that a combination of SFA initial conditions with complex weight and a trajectory model of SCTS provides the best solution for obtaining the most accurate attoclock signal [60]. The SCTS model has not been extended to the over-the-barrier ionisation (barrier-suppression regime) yet. Such an extension can be easily done as discussed above, see Eq. (2).
Recently an efficient extension and modification of the SCTS model was proposed [119]. In its orig-inal formulation the SCTS model uses the phase of the semiclassical matrix element [120][121][122] (see Refs. [123,124] for a textbook treatment), but completely disregards the pre-exponential factor of the bound-continuum transition matrix element. The influence of this pre-exponential factor was for the first time studied in Ref. [119]. The modulus of the pre-exponential factor corresponds to the mapping from initial conditions for electron trajectories to the components of the final momentum. It affects the weights of classical trajectories. The phase of the pre-exponential factor modifies the interference structures. This phase is known as a Maslov phase and can be viewed as a case of Gouy's phase anomaly, see Ref. [119]. Furthermore, a novel approach to the inverse problem applying a clustering algorithm was proposed in [119]. The modified version of the SCTS demonstrates excellent agreement with numerical solution of the TDSE for both photoelectron momentum distributions and energy spectra. It was found that the account for the preexponential factor is crucial for the quantitative agreement with the TDSE. This novel version of the SCTS can be applied not only to linearly polarised laser fields, but also to non-cylindrically-symmetric ones, e.g., bicircular laser pulses [119]. The recent semiclassical two-step model with quantum input (SCTSQI) [49] is a mixed quantumclassical approach that combines the SCTS with the numerical solution of the TDSE. To perform the synthesis of the trajectory-based approach with the TDSE, the Gabor transformation of the wave function Ψ (x, t) was used in the SCTSQI [49]. Here x 0 is the point in the vicinity of which the Gabor transform is calcu- is a Gaussian window of the width δ 0 . The quantity |G (x 0 , p x , t)| 2 describes the momentum distribution of the electron near the point x 0 at time t. This is nothing just the Husimi distribution, which can be also obtained by Gaussian smoothing of the Wigner function. In Ref. [49] the Gabor transform (4) was used in combination with the absorbing boundaries that prevent the unphysical reflections of the wave function from the grid boundary. More specifically, the Gabor transform was applied to the part of the wave function that is absorbed at every time step of the solution of the TDSE. Figure 10 shows an example of the corresponding Husimi distribution calculated at the end of a few-cycle laser pulse. This absorbed part is transformed in the ensemble of classical trajectories that is propagated using classical equations of motion. Therefore, initial positions and The absorbing boundaries must be far enough to not affect the bound part of the wave function. The SCTSQI yields quantitative agreement with quantum results [49]. What is even more important, it corrects the inaccuracies of the standard trajectorybased approaches in description of the ionisation step and circumvents the complicated problem of choosing the initial conditions. However, future work is needed to turn the SCT-SQI model in a powerful tool for studies of tunnelling. First, the model formulated for the onedimensional (1D) model atom should be generalised to the three-dimensional case (3D). To describe fine details of interference patterns accurately enough, large numbers of classical trajectories are needed in the SCTSQI. In addition to this, the ensembles of trajectories are launched at every step of the time propagation. As the result, the SCTSQI model includes all possible trajectories, and it is not always easy to distinguish between them. This hampers the understanding of the strong-field phenomena that is expected to be provided by the SCTSQI model and its future extensions. Therefore, the number of trajectories has to be reduced in the SCTSQI approach, e.g., by using more sophisticated sampling techniques. -A mask function which is absorbing the wave function over a spatial extension, such as in the SCTSQI method for example, will lead to a decreasing total probability of the wave packet. This must be mon-itored over the course of the calculation to ensure it does not introduce unwanted artefacts through the choice of position or steepness of the absorbing mask. The efficiency of this also depends on the ionisation probability which determines how much of the wave function is going to hit the absorbing boundary.

Outlook
It is evident that much remains to be done to further improve our general understanding of the tunnelling process as well as the interaction between the strong laser light and the target atom (or molecule, surface, liquid, . . . ) in order to tackle the underlying reasons for why so many approaches reach opposing conclusions. Given the lack of a clear, agreed upon definition of the onset and conclusion of tunnelling, it is perhaps unsurprising that there is also not a clear pattern between classical or quantum methods in their various predictions regarding instantaneous or finite tunnelling time, let alone numerical values. More than anything, this debate has demonstrated the need to find a common ground on which to compare the vast range of theoretical approaches and experimental setups. One of the few prevailing themes of this debate that most everyone could agree on is that a combination of classical and quantum theory is required for describing tunnelling processes in order to be able to interpret the experimental evidence.
Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: The datasets discussed in this current study are available from the corresponding author or co-authors on reasonable request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.