A probability distribution for quantum tunneling times

We propose a general expression for the probability distribution of real-valued tunneling times of a localized particle, as measured by the Salecker-Wigner-Peres quantum clock. This general expression is used to obtain the distribution of times for the scattering of a particle through a static rectangular barrier and for the tunneling decay of an initially bound state after the sudden deformation of the potential, the latter case being relevant to understand tunneling times in recent attosecond experiments involving strong field ionization.


Introduction
The search for a proper definition of quantum tunneling times for massive particles, having well-behaved properties for a wide range of parameters, has remained an important and open theoretical problem since, essentially, the inception of quantum mechanics (see, e.g., [1,2] and references therein). However, such tunneling times were beyond the experimental reach until recent advances in ultrafast physics have made possible measurements of time in the attosecond scale, opening up the experimental possibility of measuring electronic tunneling times through a classically forbidden region [3][4][5][6] and reigniting the discussion of tunneling times. Still, the intrinsic experimental difficulties associated with both the measurements and the interpretation of the results have, so far, prevented an elucidation of the problem and, in fact, contradictory results persist, with some experiments obtaining a finite nonzero result [3,6] and others compatible with instantaneous tunneling [4]. It should be noticed that the similarity between Schrödinger and Helmholtz equations allows for analogies between quantum tunneling of massive particles and photons [7], and a noninstantaneous tunneling time is supported by this analogy and experiments measuring photonic tunneling times [8], as well as by many theoretical calculations based on both the Schrödinger (for reviews see, e.g., [1,2]) and the Dirac equations (e.g., [9][10][11][12][13][14][15][16]).
The conceptual difficulty in obtaining an unambiguous and well-defined tunneling time is associated with the impossibility of obtaining a self-adjoint time operator in quantum mechanics [17], therefore leading to the need for operational definitions of time. Several such definitions exist, such as phase time [18], dwell time [19], the Larmor times [20][21][22][23], and the Salecker-Wigner-Peres (SWP) time [17,24], and in some situations these lead to different, or even contradictory, results. This is not surprising, since by their own nature operational definitions can only describe limited aspects of the phenomena of tunneling, and it is unlikely that any one definition will be able to provide a unified description of the quantum tunneling times in a broad range of situations. Nevertheless, it remains an important task to obtain a welldefined and real time scale that accurately describes the recent experiments [3][4][5][6][25][26][27].
It is important to notice that the time-independent approach to tunneling times (i.e., for incident particles with sharply defined energy), which comprises the vast majority of the literature, is ill-suited to accomplish the above-mentioned goal, since it ignores the essential role of localizability in defining a time scale [23,28]; see, however, [29], which applies the time defined in [30] to investigate the half-life of -decaying nuclei. A few works (e.g., [23,28,31,32]) address the issue of localizability and, consequently, arrive at a probabilistic definition of tunneling times (that is, an 2 Advances in High Energy Physics average time). In particular, in [28] the SWP clock was used to obtain an average tunneling time of transmission (reflection) for an incident wave packet, and such time was employed to investigate the Hartman effect [33] for a particle scattered off a square barrier and it was shown that it does not saturate in the opaque regime [28,34].
The tunneling time scales considered in [23,28,31] involve taking an average over the spectral components of the transmitted wave packet and, thus, obscure the interpretation of the resulting average time. In this paper, we take as a starting point the real-valued average tunneling time obtained in [28], using the SWP quantum clock, and obtain a probability distribution of transmission times, by using a standard transformation between random variables. In addition to providing a more accurate time characterization of the tunneling process, this should provide a clearer connection with the experiments (which measure a distribution of tunneling times; see, e.g., Figure 4 in [3]). It is worth noting that some approaches using Feynman's path integrals address the problem of obtaining a probabilistic distribution of the tunneling times (see, e.g., [35]). However, these methods in general result in a complex time (or, equivalently, multiple time scales), and some arbitrary procedure is needed to select the physically meaningful real time a posteriori.
After obtaining a general formula for the distribution of tunneling times, which is the main result of this work, we apply it to two specific cases. First, to illustrate the formalism in a simple scenario, we consider the situation of a particle tunneling through a rectangular barrier. Then, we consider a slight modification of the model proposed in [36] for the tunneling decay of an initially bound state, after the sudden deformation of the binding potential by the application of a strong external field; the modification considered here allows us to investigate the whole range of possibilities for the tunneling times, without having an "upper cutoff ", as is the case in the original model. Finally, some additional comments on the results are reserved for the last section.

The SWP Clock's Average Tunneling Time
We start by briefly reviewing the time-dependent application of the SWP clock to the scattering of a massive particle off a localized static potential barrier in one dimension (for details see [28]) which is appropriate, since it follows from the threedimensional Schrödinger equation for this problem that the dynamics is essentially one-dimensional [3].
The SWP clock is a quantum rotor weakly coupled to the tunneling particle and that runs only when the particle is within the region in which ( ) ̸ = 0, where ( ) is the potential energy. The Hamiltonian of the particle-clock system is given by (we use ℏ = 2 = 1, where is the particle's mass) [17] = − where P( ) = 1 if ( ) ̸ = 0 and zero otherwise. The clock's Hamiltonian is = − ( / ), where the angle ∈ [0, 2 ) is the clock's coordinate and = 2 /(2 + 1) is the clock's angular frequency, with being a nonnegative integer or halfinteger giving the clock's total angular momentum, and is the clock's resolution. The weak coupling condition amounts to assume that is large, in such a way that the clock's energy eigenvalues, ≡ (− < < ), are very small compared to the barrier height and the particle's energy. It is assumed that, at = 0, well before it reaches the barrier, the particle is well-localized far to the left of the barrier and the wave function of the system is a product state of the form where ( ) is the particle's initial state, represented by a wave packet centered around an energy 0 , and the clock initial state is assumed to be "in the zero-th hour" [17] where ( ) = e / √ 2 are the clock's eigenfunctions corresponding to the energy eigenvalues .
The state V 0 ( ) is strongly peaked at = 0, thus allowing the interpretation of the angle as the clock's hand, since for a freely running clock the peak evolves to , where is the time measured by the clock [17]. Since here clock and particle are coupled according to (1), when the particle passes through the region ( ) ̸ = 0 it becomes entangled with the clock, with the wave function for the entire system given by where is the incident particle's energy, = √ , and ( ) is the Fourier spectral decomposition of the initial wave packet ( ) in terms of the free particle eigenfunctions (we are assuming delta-normalized eigenfunctions). The functions ( ) ( ) satisfy a time-independent Schrödinger equation with a constant potential in the barrier region. Outside the potential barrier region and for a particle incident from the left, the (unnormalized) solution ( ) ( ) of the timeindependent Schrödinger equation is given by [28] where ( ) ( ) [ ( ) ( )] stands for the transmission (reflection) coefficient, and it is assumed, without loss of generality, that the potential is located in the region − < < . Considering only the transmitted solution in (5) and substituting it into the time-dependent solution (4), it can be shown that for weak coupling Advances in High Energy Physics is the stationary transmission clock time corresponding to the wave number component [17,37]. The transmission coefficient ( ) corresponds to the stationary problem in the absence of the clock. For tunneling times one is interested only in the clock's reading for the postselected asymptotically transmitted wave packet. Thus, tracing out the particle's degrees of freedom, the expectation value of the clock's measurement can be defined, resulting in the average tunneling time [28] where = 1/ ∫ | ( ) ( )| 2 is a normalization constant and ( ) is the probability density of finding the component in the transmitted wave packet. Similar expressions can be obtained for the reflection time.

The Tunneling Times Distribution
An important aspect of the average tunneling time considered in the previous section is that it emphasizes the probabilistic nature of the tunneling process. However, since the average in (8) is over the time taken by the spectral components of the wave packet, it does not lend itself to an easy interpretation, given the spectral components of the wave packet tunnel with different times. Thus, instead of (8), one would rather obtain an average over (real) times of the form where ( ) stands for the probability density for observing a particular tunneling time for the asymptotically transmitted wave packet. This can easily be achieved by noticing that in probability theory (8) and (9), which must be equal, are related by a standard transformation between the two random variables and through a function ( ). It follows that the probability distribution of times is given by which, in essence, is the statement that all the -components in the transmitted packet for which ( ) = must contribute to the value of ( ) with a weight ( ). Finally, using the properties of the Dirac delta function (specifically, we use the fact that ( ( )) = ∑ ( ( − )/| ( )|), where { } is the set of zeros of the function ( ) and the prime indicates a derivative with respect to the independent variable), we obtain where { ( )} is the set of zeros of the function ( ) ≡ ( )− and is the derivative of ( ) with respect to . A similar definition of the distribution of tunneling times given in (10)-(11) can be obtained for any time scale which is probabilistic in nature, that is, of the form (8).
Although several other probabilistic tunneling times exist in the literature (e.g., [23,31,32,35]), the SWP clock has proven to yield well-behaved real times both in the time-independent [17,37,38] and time-dependent approaches [28,34,39] and it provides a simple procedure to derive the probabilistic expression (8). In addition, the role exerted by circularly polarized light in attoclock experiments [3,25] seems to provide a natural possibility for interpretation in terms of the SWP clock.
As will be illustrated below, for the simple application of this formalism to the problem of a wave packet scattered off a rectangular potential barrier, the distribution of times (10)-(11) cannot, in general, be obtained analytically even for the simplest cases, except in trivial cases such as for a single Dirac delta potential barrier [40][41][42], in which case ( ) = 0 and ( ) = ( ) ∫ ( ). It should also be noticed that, despite the fact that the derivation of the previous section leading to (8) and, thus (10)-(11), assumed a scattering situation, these expressions can be shown to be valid for any situation involving preselection of an initial state localized to the left of a potential "barrier" followed by postselection of an asymptotic transmitted wave packet. This allows us to obtain the distribution of times for a model that simulates the tunneling decay of an initially bound particle by ionization induced by the sudden application of a strong external field; the model considered below is a variant of that introduced in [36].

The Distribution of Tunneling Times for a Rectangular Barrier
As a first illustration of the formalism developed above, let us consider a rectangular barrier of height 0 located in the region ∈ (− , ). The particle's initial state 0 ( ) ≡ ( , = 0) is assumed to be a Gaussian wave packet where the parameters 0 , , and 0 are chosen such that the wave packet is sharply peaked in a tunneling wave number 0 = √ 0 < √ 0 and is initially well-localized around = 0 , far to the left of the barrier; in the calculations that follow we take 0 = −8 , such that at = 0 the probability of finding the particle within or to the right of the barrier is negligible. The transmission coefficient ( ) and the spectral function ( ) are well-known and given by Advances in High Energy Physics where = √ 0 − 2 . The stationary transmission clock time (7) is [22,28] with tunneling times corresponding to real values of (i.e., 0 > 2 ). Figure 1 shows a plot for the stationary transmission times ( ), the distribution of wave numbers ( ) in the transmitted wave packet, and the distribution | ( )| 2 of wave numbers (momenta) in the incident packet, for two values of the barrier width. For the chosen parameters and barrier widths both the incident and the transmitted wave packets have an energy distribution very strongly peaked in a tunneling component (in the bottom plot of Figure 1 the barrier is much more opaque than that in the top plot and we can observe that-despite being with a negligible probability for the parameters chosen for this plot-in this situation some above-the-barrier components start to appear in the distribution of the transmitted wave packet. So, in order to consider mainly transmission by tunneling we must restrict the barrier widths to not too large ones). We also observe the very well-known fact that the transmitted wave packet "speeds up" when compared to the incident particle [28]. As a general rule, the larger is the barrier width (i.e., the more opaque is the barrier), the greater is the translation of the central component towards higher momenta. In what concerns the off-resonance stationary transmission time, it initially grows with the barrier width, and saturates for very opaque barriers (the Hartman effect); on the other hand, it presents peaks at resonant wave numbers that grow and narrow with the barrier width; for a detailed discussion see [28]). Figure 2 shows plots of the probability distribution ( ) of the tunneling times according to (10)-(11), corresponding to both the barrier widths shown in Figure 1 [to obtain these plots we used a Monte Carlo procedure to generate a large number of outcomes from the distribution ( ), which afterwards were transformed into values by using the function = ( )]. The vertical grey lines in these plots correspond to the time the light takes to cross the barrier distance. It is observed that for the two distributions shown in Figure 2 the probability to observe superluminal tunneling times is negligible. It is also observed that these distributions have a shape that resembles that of the distribution, albeit with a more pronounced skewness. This shape could be inferred from Figure 1 and from (11), since T ( ) grows very smoothly in the region were ( ) is nonvanishing. Furthermore, a comparison between the two plots in Figure 2 shows that the tunneling times do not grow linearly with the barrier width and, therefore, the distribution in the bottom plot of Figure 2 is "closer" to the light time than the distribution shown in the top plot; [28] already observed that for intermediate values of barrier widths the average transmission time-corresponding to the mean of the distribution -reaches a plateau.

Distribution of Ionization Tunneling Times
In this section we obtain a distribution for tunneling times for a particle that is initially in a bound state of a given binding potential. The potential is then suddenly deformed in such a way that the particle can escape from the initially confining region by tunneling. The model considered here is a slight modification of that proposed by Ban et al. [36] to simulate, Advances in High Energy Physics in a simple scenario, key features of the decay of a localized state by tunneling ionization induced by the application of a strong external field with a finite duration.
In [36], for < 0, the particle is in an eigenstate of a semiinfinite square-well potential 1 ( ), and, therefore, it cannot decay by tunneling. At = 0 the potential is suddenly deformed to 2 ( ), such that the particle can now tunnel through the potential barrier; it is assumed that the wave function does not change during the sudden change of the potential. Finally, after a finite time 0 the potential returns to its original configuration, 1 ( ), and tunneling terminates. The cutoff time 0 mimics the natural upper bound for tunneling times measured in recent attoclock experiments (see, e.g., [3,6] and references therein), since the opening and closing of the tunneling channel in these experiments occur in intervals of half the laser field's period.
Here, we deviate from [36] by setting 0 → ∞; i.e., once deformed the potential does not return to its original form and, after a long enough time, the particle will be transmitted with unit probability; thus, by eliminating the cutoff (which is just an experimental limitation) we are able to explore the whole range of possibilities for the ionization tunneling time. In addition, for ≥ 0, the particle is assumed to be coupled to a SWP quantum clock running only in the region ( , ), so that the clock's readings for the asymptotic transmitted wave packet give the time the particle spent within the barrier after = 0. Following [36], we assume that for < 0 the particle is in the ground state of the potential 1 ( ), whose stationary wave function is given by where is a normalization constant, 0 = √ 0 , 0 is the ground state energy, and 0 = √ 0 − 2 0 . It is also assumed, as in [36], that immediately after the sudden deformation of the potential from 1 ( ) to 2 ( ), at = 0, the wave function does not change. However, for ≥ 0 the particle state, which is no longer an energy eigenstate, is given by a superposition of the energy eigenstates ( ) ( = √ ) of the potential 2 ( ), i.e., [36] where with where = √ 0 − 2 and the coefficients ( ), ( ), ( ) and the phase Ω( ) are determined by the usual boundary conditions at = and = and are such that the normalization ⟨ ( ), ( )⟩ = ( − ) holds [36]. From the above expressions it follows that, without any loss of generality, we can take ( ) and all the eigenfunctions (21) to be real.

Advances in High Energy Physics
In order to consider the coupling with the SWP clock for times ≥ 0 we proceed as follows. At = 0 the system particle+clock is described by the product state ( , 0)V 0 ( ), where ( , 0) is the state (19) and V 0 ( ) is the initial clock state given by (3). After = 0 the particle and the clock states become entangled. For the procedure of postselection of the asymptotically transmitted wave function we notice that the role of the transmission coefficient for the wave function (21) is played by √2/ e (− +Ω ( ) ( )) , where the superscript indicates the weak coupling with the clock. The right moving asymptotic wave packet representing the coupled system formed by the transmitted particle and the clock is where, as before, ( ) = −( Ω ( ) ( )/ ) =0 = −(1/ 2 )( Ω/ ) [with quantities without the subscript "( )" representing the limit → 0]. By following the same steps described in [28], we trace out the clock's degree of freedom in the asymptotic transmitted wave packet in order to obtain the distribution ( ) of the wave numbers for the asymptotically transmitted wave packet, which in this case is simply given by i.e., the probability to find a wave number in the asymptotic transmitted wave packet is the same as in the initial state, which is as expected, since after a long enough time the initial wave packet will be transmitted with probability unit, as mentioned earlier.
The general behavior of ( ) and ( ) is illustrated in Figures 3 and 4, corresponding to two barriers with different opacities ( − = 2 and 4, respectively). These plots show, as expected, that the distribution ( ) is strongly peaked at the wave number 0 , corresponding to the energy of the initially bound state and is negligible for nontunneling components. For tunneling wave numbers ( < √ 0 ) the function ( ) is also strongly peaked at the same wave number 0 , which corresponds to a local maximum (for nontunneling wave numbers there are several other resonance peaks). From (11) we would expect that the peaks in the tunneling times distribution ( ) would occur for times = ( ) corresponding to values of for which ( ) ≈ 0-which occur at points of local maxima and minima of the function ( )-and corresponding to nonnegligible ( ). Therefore, from the plots in Figures 3 and 4 one could expect the first peak of the tunneling time distribution ( ) at ≈ 0.105 . . (the local minimum of ( ), which is similar for both barrier widths, since nonresonant times ( ) change little with the barrier width for opaque barriers, as is the case in Figures 3 and 4); a second peak in ( ) is expected to occur around the local maximum of ( ), which corresponds to ≈ ( 0 ) (this local maximum-corresponding to resonant wave numbers-changes significantly with the barrier widths; see, e.g., [28]). On the other hand, peaks in with the initial state given by (18). Bottom: close view of the above plot for small times. The vertical grey lines in the plots correspond to = 0 and = √ 0 . The regions in which ( ) ≈ 0 (around the local maximum and minimum of ( )) correspond to times ≈ ( 0 ) and ≈ 0.105 . . Rydberg atomic units were used in all the plots.
( ) coming from local maxima (resonances) and minima associated with nontunneling values of are suppressed, since ( ) ≈ 0 in these cases. Figure 5 confirm these claims. For both barrier widths considered, the distribution of tunneling times is "U" shaped, having peaks at the times corresponding to the local maxima and minima of the stationary time ( ) inside the tunneling region. It should be observed that the larger is the barrier width, the broader is the tunneling time distribution, due to the strong increase of the resonant tunneling time with the barrier width. region, which in both plots corresponds to almost the same value ≈ 0.105 . ≈ 5.1 attoseconds. The top plot of Figure 6 shows that for the less opaque barrier there exists a (very small) probability to observe a superluminal tunneling time. Even if this possibility cannot be precluded in principle (see, e.g., [16]), in the present case the possibility of emergence of such small times was expected, since at = 0 there was a significant portion of the wave packet (roughly 27%) penetrating the whole distance of the barrier, and this has an important contribution to the emergence of small times in the clock's readings associated with the transmitted particle. On the other hand, the top plot of Figure 7 shows that for the thicker barrier the probability for superluminal times is negligible; the portion of the wave packet already inside the barrier at = 0 is the same (∼27%), but the wave packet penetrates proportionally a smaller distance inside the barrier and, thus, it does not contribute in a significant way to the  Figure 5: Distributions of tunneling decay times ( ) through the barrier of the potential 2 ( ) for the initial bound state 0 ( ) given by (18). The histograms were built by using the Monte Carlo procedure described in Figure 2 and in the main text. The vertical grey lines indicate percentiles of the distribution (the first and the last correspond to 1% and 99%, the remaining ones range from 5% to 95%, in steps of 5%); the three thick vertical lines indicate the first quartile (percentile 25%), the median, and the third quartile emergence of very small times in the clock readings. We note that the introduction of the cutoff 0 , as in [36], would result in a time distribution similar to the truncated distributions shown in the top plots of Figures 6 and 7. It is also worth observing that, for small times, the distributions obtained here resemble qualitatively those in Figure 4 of [3], except for the presence of several peaks at discrete values of the time in the latter. The considerations above, relating the peaks of the distribution of clock times ( ) to the local maxima and minima of the stationary time ( ) and the magnitude of distribution ( ) in the neighborhood of these points, suggest a scenario in which such multiple peaks at discrete values of time can appear in the distribution ( ) of transmission times. Indeed, if abovethe-barrier wave numbers had a significant contribution to the initial wave packet, then the several local maxima and minima present in the vicinities of the resonant nontunneling components will also contribute in a significant way to build multiple peaks in the distribution of transmission times; these peaks, however, could not be associated with the tunneling process. We can consider such a scenario by choosing as the initial state a tightly localized state given by ( , 0) = √ 2 sin 0 , with 0 = and the barrier parameters = 1, = 2, and 0 = 11, in Rydberg atomic units. In this situation the initial wave function is perfectly confined to the left of the barrier (0 < < 1), and above-the-barrier components contribute in a significant way to build the wave packet, as can be seen from ( ) in the top plot of Figure 8 (in this case the probability of finding a nontunneling component in the wave packet is approximately 75%). In this plot we can also observe that all the local maxima and minima of ( ) shown occur in neighborhoods of wave numbers for which ( ) is nonnegligible; therefore, all these local maxima and minima contribute significantly to build multiple peaks in the distribution of transmission times ( ). The middle and the bottom plots of Figure  peaks in the bottom plot are associated with nontunneling components).

Conclusions
Taking as a starting point the probabilistic (average) tunneling time obtained in [28] with the use of a SWP clock [17,24,37], we obtained a probability distribution of times (10)- (11). An important advantage of using the SWP clock, in addition to those already mentioned, is that by running only when the particle is inside the barrier it allows us to address the concept of tunneling time in a proper way, since the time spent by the particle standing in the well before penetrating the barrier is not computed. A clear advantage of having a probability distribution of transmission (tunneling) times is that, in addition to the usual expectation value, we can obtain all the statistical properties of this time, such as its most probable values (peaks of the distribution), the dispersion around the mean value, and the probability to observe extreme outcomes (superluminal times, for instance).
As an initial test, the distribution of times (10)-(11) was applied to the simple problem of a particle tunneling through Advances in High Energy Physics a rectangular barrier. Unsurprisingly, it revealed behavior similar to that already known from previous works using a distribution of wave numbers (momentum)-see, e.g., [28]-although, using ( ) these conclusions are much more transparent. For example, one could answer the question about the possibility of superluminal tunneling by direct calculation from the probability distribution ( ). In the nonstationary case-which is the correct to address this question-this problem is usually answered by considering just the average tunneling time. But, given its probabilistic nature, an answer based only on the average time may not be satisfactory, especially if the dispersion of the distribution of tunneling times is large, which is often the case when one deals with well-localized particles, as suggested by the two situations addressed in the present work.
As a main application of (10)-(11), we considered a slight modification of the problem considered in [36] to model strong field ionization by tunneling. The modification considered here was the elimination of the cutoff time that was introduced in [36] to simulate the upper bound that arises in attoclock experiments [3,6] due to the opening and closing of the tunneling channel, naturally associated with the oscillations in the laser field intensity. This cutoff is not a fundamental requirement, but rather it is associated with the experimental methods employed-in any case, its implementation is rather trivial, since it just truncates the distribution of times. The consideration of the full range of the distribution of times allowed us to show that an important contribution to ( ) comes from very large times associated with the resonance peaks in the tunneling region; these very long tunneling times occur with a probability comparable to very short ones, thus having an important impact on the average tunneling times and, therefore, cause difficulties when comparing theoretical predictions based on an average time with the outcomes of experiments presenting a natural cutoff in the possible time measurements. In particular, in the attoclock experiments the relevant measure is often associated with the peak of the tunneling time, which may be promptly identified once one knows the probability distribution for all possible times. A remark is in place; the distribution of times proposed here, built on the SWP clock readings, refers to the time the particle dwells within the barrier, while the tunneling times often measured in recent attoclock experiments actually refer to exit times [43].
In sum, the approach introduced above and resulting in (10)-(11) builds upon the already (conceptually) well tested SWP clock to provide a real-valued distribution of times that, in the simple models considered here, was demonstrated to have physically sound properties and, in fact, (rough) similarities with the time distribution obtained in recent experiments [3], therefore warrantying further investigation with more realistic potentials.

Data Availability
The numerical data used in the construction of the histograms in this article were generated by a Monte Carlo procedure, whose details were described in the main text.