Analysis of wavepacket tunneling with the method of Laplace transformation
Natascha Riahi
quant-ph

Polarizability Extraction for Waveguide-Fed Metasurfaces
    David Smith
    July 14, 2017
========================================================


 We use the method of Laplace transformation to determine the dynamics of a wave
packet that passes a barrier by tunneling. We investigate the transmitted wave packet
and find that it can be resolved into a sequence of subsequent wave packages. This result
sheds new light on the Hartman effect for the tunneling time and gives a possible explanation for an 
experimental result obtained by Spielmann et. al.




§ INTRODUCTION


There are several definitions
 of tunneling times (<cit.>,<cit.>,<cit.>,<cit.>,<cit.>,<cit.>,<cit.>,<cit.>)
 and the discussion about their meaning is still 
ongoing
(<cit.>,<cit.>,<cit.>,<cit.>). Recently, 
experiments with atoms stimulated by ultrashort, 
infrared laser pulses, the so-called attoclock experiments,
reinforced the interest in the prediction of tunneling times (<cit.>,<cit.>,<cit.>).

In this article we  investigate the dynamics of a wavepacket that tunnels through a barrier.
The tunneling time we determine is the
 so-called group delay or phase time (<cit.>,<cit.>). It can be characterized as the time interval
 between the moment the peak of a freely evolving incoming  wave packet would reach the barrier and the arrival of the peak of the transmitted
 wave packet at the end of the barrier.
 One of the oldest results for this tunneling time was obtained by Mac Coll 
(<cit.>) in 1932 by the application of the stationary phase method. He concluded
that there was ”no appreciable delay in the transmission of the packet through the barrier”. These calculations were later 
refined by Hartmann (<cit.>,<cit.>), 
who found a finite delay time. For thicker barriers this delay
time becomes independent of the thickness of the barrier and tends to a fixed value
which is known as the Hartmann effect. This saturation  was also found
for other definitions of tunneling time (<cit.>), whereas in the framework of fractional quantum mechanics a decreasing of the 
tunneling time with barrier width was obtained(<cit.>).  A formal analogy between the Schrödinger and the 
Helmholtz equation (see for instance <cit.>)made it possible to test predictions for the tunneling time with optical experiments. 
The results confirmed the saturation of the tunneling time (<cit.>,<cit.>, <cit.>,see <cit.> for more references ) 
though a more detailed discrepancy between experiment and calculations remained open in <cit.>.

 Since the Hartmann effect seems to threaten Einstein causality, because superluminal velocities can be infered from transmission times saturating with barrier width ,  an objection was made that only 
 the modes above the barrier energy may contribute to the transmitted wave packet
 (<cit.>, <cit.>). But it was argued in (<cit.>)with reference to numerical calculations and experimental results that this effect can not be used as an explanation of the barrier tunneling phenomenon  (see also <cit.> for a more detailed discussion and further references).

The results of the attoclock experiments with strong laser fields that lower the Coulomb potential  can be best explained by using a 
tunneling time probability amplitude constructed with Feynnman path integrals (<cit.>, <cit.>), whereas several other  definitions of tunneling times do 
 not agree with the data.
 
 At present there exists no unified  formalism for the calculation of tunneling times for different experimental settings. 
 In the case of wavepacket tunneling through static barriers the phase time, keeping track of the peak is the best candidate for the tunneling time, 
 and in good agreement with the experimental results(<cit.>,<cit.>, <cit.>).


 
The phase time is in general obtained by representing the solution of the 
 time-dependent Schrödinger equation as the integral over stationary solutions for different energies. 
Under the assumption that the
energy distribution of the initial wave packet is sufficiently peaked the application of the 
stationary phase
method yields a time for the arrival of the peak of the wavepacket at the end of the barrier. 

In this article we use the method of Laplace transformation
<cit.> to determine the wavepacket dynamics.Instead of immediately applying an approximation for oscillatory integrals, we first obtain exact solutions for the wavefunction at each side and inside  the barrier, which we then simplify assuming that the initial wavepacket is sufficiently 
 peaked in the momentum 
representation. Since we are free to choose explicitely the initial wavepacket
in position space (which we do not have to construct as superposition of eigenstates), it is
technically much easier to start with a wavepacket located exclusively at the left side of the barrier. 
Furthermore we give estimates for the consequences of the applied approximations and derive a consistency condition in the case of a special class of initial wavepackets.

We find that what comes out of the barrier is not a 
single
wave packet, but infinitely many, one after the other. The time each of these  wavepackets needs 
for
the tunneling 
is completely independent of the thickness of the barrier. But for thicker barriers, the later
wavepackets are much more attenuated and so there remains a single transit time that equals the result Hartmann obtained as 
the upper limit for the tunneling time. However in our approach the tunneling decreases before it reaches the limiting value
for thick barriers. This is in accordance with the results of the optical experiments performed by
Spielmann et al. <cit.>.
 Bernardini (<cit.>) investigated the transmission of wavepackets with energies above the height of the barrier.
He also finds that there is not only one reflected and one transmitted wave
packet, but many of them that leave the barrier. This result was obtained by 
the decomposition
of the transmission and reflection coefficients of the stationary solution in  
 infinite series.  Our result is the counterpart  in the tunneling regime. 



§ THE SOLUTION OF THE TIME-DEPENDENT SCHRÖDINGER EQUATION


A finite barrier 
 is described by the potential



    [    V(x)=0               x ≤ 0;    V(x)=V           0 < x < d;    V(x)=0           x ≥ d   . ]
 


The dynamics of a wave packet ψ(x,0)  is determined by the Schrödinger equation 


    -ħ^2/2 m∂^2 ψ(x,t)/∂ x^2
     = i  ħ ∂ψ(x,t)/∂ t            x < 0  ,
        -ħ^2/2 m∂^2 ψ(x,t)/∂ x^2+V ψ(x,t)
     = i  ħ ∂ψ(x,t)/∂ t             0 < x < d  ,
       -ħ^2/2 m∂^2 ψ(x,t)/∂ x^2
     = i  ħ ∂ψ(x,t)/∂ t            x > d  ,

 
 where ψ(x,t) is supposed to be continuously differentiable everywhere and square integrable .
 

We will further assume that the initial wave packets are located
at the left side of the barrier

 
    ψ(x,0)=0          x ≥ 0 .

 

 

We apply the method of Laplace transformation, already introduced in (<cit.>). 

The Laplace transformed wavepacket φ(x,s)


    φ(x,s) = ℒ(ψ(x,t)) = ∫_0^∞ψ(x,t)  e^-s t dt    .



obeys the transformed equations

    -ħ^2/2 m∂^2 φ(x,s)/∂ x^2 φ(x,s) = i ħ  s  φ(x,s) - i ħ ψ(x,0)            x < 0  ,
       -ħ^2/2 m∂^2 φ(x,s)/∂ x^2 + V  φ(x,s) = i ħ  s  φ(x,s) - i ħ ψ(x,0)              0 < x < d  ,
       -ħ^2/2 m∂^2 φ(x,s)/∂ x^2 φ(x,s) = i ħ  s  φ(x,s) - i ħ ψ(x,0)             x > d  .

 
 The solution of (<ref>) is determined by the method of variation of constants:
 





    φ(x,s)   =√(m/2 i s ħ){u_1(x,s)∫_-d^xu_2(y,s) 
    ψ(y,0) dy    -  .   .
     u_2(x,s) 
    ∫_-d^xu_1(y,s) ψ(y,0) dy 
    }
          
          + α(s) u_1(x,s)   + β(s) u_2(x,s) 
           x < 0     
          
    φ(x,s)   = γ(s) u_3(x,s)   + δ(s) u_4(x,s)
           0 <  x < d     
    φ(x,s)   = μ(s) u_1(x,s)   + ν(s) u_2(x,s)
           x > d.



The functions u_1(x,s),u_2(x,s),u_3(x,s),u_4(x,s) are the solutions of the homogeneous
equations corresponding to (<ref>)



    u_1(x,s)=e^i√(2 m s i/ħ) x      u_2(x,s)=e^-i √(2 m s i/ħ) x      0 <  x < d
    
        u_4(x,s)=e^-i √(2 m s i/ħ-2 m V/ħ^2) x     u_3(x,s)=e^i√(2 m s i/ħ-2 m V/ħ^2) x      x <  0,x > d.



Since φ(x,s) must vanish for x→±∞, we find

    α(s)=0 ,  ν(s)=0 .


If we evaluate φ(x,s) and its first derivative at x=0 and x=d, we obtain, imposing continuous differentiability


    β(s)=√(m/2 i ħ s)(I_1-I_2)
    +√(m/2 i ħ s)·2  I_2/1+√(1-V/i ħ s)·ρ(s)-e^-2 i d √(2 m s i/ħ-2 m V/ħ^2)/ρ^2(s)-e^-2 i d √(2 m s i/ħ-2 m V/ħ^2)
       γ(s)=
    -√(m/2 i ħ s)·2  I_2/1+√(1-V/i ħ s)·e^-2 i d √(2 m s i/ħ-2 m V/ħ^2)/ρ^2(s)-e^-2 i d √(2 m s i/ħ-2 m V/ħ^2)
       δ(s)=√(m/2 i ħ s)·2  I_2 ρ(s)/1+√(1-V/i ħ s)·1/ρ^2(s)-e^-2 i d √(2 m s i/ħ-2 m V/ħ^2) ,
       μ(s)=√(m/2 i ħ s)·2  I_2/1+√(1-V/i ħ s)
    e^- i d √(2 m s i/ħ)·1-ρ(s)/ρ^2(s)-e^-2 i d √(2 m s i/ħ-2 m V/ħ^2)
    e^-2 i d √(2 m s i/ħ-2 m V/ħ^2) ,

where we have introduced the abbreviations

    I_1=∫_-∞^0e^i √(2 m s i/ħ)yψ(y,0)dy ,
      
    I_2=∫_-∞^0e^-i √(2 m s i/ħ)yψ(y,0)dy

and also used

    ρ(s)=2/ 1+√(1-V/ħ s i)-1 .


Inserting this result into (<ref>) and applying the  series expansion 

    1/1-ρ^2(s)e^2 i d √(2 m s i/ħ-2 m V/ħ^2)=∑_k=0^∞
     (ρ(s))^2 ke^2 i d k√(2 m s i/ħ-2 m V/ħ^2)

yields 



    φ(x,s)=
     √(m/2 s i ħ){∫_-∞^0
    e^i√(2 m s i/ħ) |x-y|ψ(y,0) dy+ 
    ∫_-∞^0
    e^-i√(2 m s i/ħ) (x+y)ψ(y,0) dy
    ·ρ(s)} + 
           √(m/2 s i ħ)∑_l=0^∞∫_-∞^0
    e^2 i d (l+1)√(2 m s i/ħ-2 m V/ħ^2)
    -i (x+y)√(2 m si/ħ)ψ(y,0)dy
    · (ρ^2(s)-1)ρ^2l+1(s)
    
              x < 0 
       φ(x,s)=
     √(m/2 s i ħ)∑_l=0^∞∫_-∞^0
    e^ i  (2 d l+x)√(2 m s i/ħ-2 m V/ħ^2)
    -i y √(2 m si/ħ)ψ(y,0)dy
    · (ρ(s)+1)ρ^2l(s) -
           √(m/2 s i ħ)∑_l=0^∞∫_-∞^0
    e^ i  (2 d (l+1)-x)√(2 m s i/ħ-2 m V/ħ^2)
    -i y√(2 m si/ħ)ψ(y,0)dy
    · (ρ(s)+1)ρ^2l+1(s) 
              0 ≤ x ≤ d  ,
       φ(x,s)=
     √(m/2 s i ħ)∑_l=0^∞∫_-∞^0
    e^i d (2l+1)√(2 m s i/ħ-2 m V/ħ^2)
    +i (x-d-y)√(2 m si/ħ)ψ(y,0)dy
    · (-ρ^2(s)+1)ρ^2l(s) 
              x > d  .




We introduce the shifted momentum representation 

    f(Q,p)=1/√(2 πħ)∫_-∞^∞ψ(x+Q,0) e^-i p x/ħ dx = 
     e^i p Q/ħf(p) ,

where 

    f(0,p) = f(p)

denotes the representation of the wave function in momentum space.

Applying the abbreviations


    a_l(t)=ℒ^-1{ρ^2 l+1(s)(1-ρ^2(s))} ,   
    b_l(t)=ℒ^-1{ρ^2 l(s)(1+ρ(s))}
        c_l(t)=ℒ^-1{ρ^2 l+1(s)(1+ρ(s))} ,   
    g_l(t)=ℒ^-1{ρ^2 l(s)(1-ρ^2(s))} ,

we find proceeding as for the asymmetric square well in <cit.> 
 for the inverse Laplace transform of (<ref>)



    ψ(x,t)=
    κ/√(2 π t  i)∫_-∞^0  e^i (x-y)^2 κ ^2/2 tψ(y,0) dy
    +
    ∫_0^tκ/√(2 π (t-τ)  i)∫_-∞^0e^i (x+y)^2 κ ^2/2 (t-τ)ψ(y,0) r(τ)dy  dτ + 
         ∑_l=0^∞∫_0^t∫_-∞^∞
    K(2d( l+1),p,t-τ)
    f(-x,p)dp · a_l(τ)dτ
            x < 0  , 
       ψ(x,t)=
    ∑_l=0^∞∫_0^t∫_-∞^∞
    K( 2 d l+x),p,t-τ)
    f(p)dp · b_l(τ)dτ + 
         ∑_l=0^∞∫_0^t∫_-∞^∞
    K(d (2l+1)-x),p,t-τ)
    f(p)dp · c_l(τ)dτ 
            0 < x < d  , 
       ψ(x,t)=
    ∑_l=0^∞∫_0^t∫_-∞^∞
    K(d(2 l+1),p,t-τ)
    f(x-d,p)dp · g_l(τ)dτ 
            x > d  ,


where K(x,p,t) is defined by

    K(x,p,t)= 1/2 √(2 πħ) e^- i V t/ħ-i t q^2/2 ħ m  
         {
     e^-i x q/ħ[-i √(2 m i/ħ  t)x/2-i √(i  t/2 ħ m) q]
     +e^i x q/ħ[-i √(2 m i/ħ  t)x/2+i √(i  t/2 ħ m) q]
    } ,
 
 and we used κ= √(m/ħ) and q=√(p^2-2 m V).
 


§ TUNNELING OF WAVE PACKETS




In order to investigate the tunneling process, we
 consider  a wave packet that is represented in momentum space by 
 


    f(p)=e^-i p x_0/ħ F(p-p_0)   ,


where F(p) fulfills

    ∫_-∞^∞ F^*(p) p  F(p) dp=0  ,    ∫_-∞^∞  F^*(p)F'(p) dp=0 .

.
The expectation values of position and momentum are then given by

    ⟨x̂⟩=x_0 ,
        ⟨p̂⟩=p_0 .

We  assume that the wavepacket is  concentrated within a region 0 < p <p_max <√(2 m V).
so that we can use the approximation

    f(p) ≈ 0      p > p_max      p < p_min .

Moreover it should be sufficiently peaked around the momentum expectation value to justify the approximation

    F(p-p_0) (p-p_0)
     ≈ 0 .

Finally the difference √(2 m V)-p_m should be big enough
to ensure 

    1/√(1-p^2/2 m V) = O(1)     p_min ≤  p ≤ p_max .

In Appendix <ref> it is shown that these assumptions about the momentum distribution are compatible with the requirement (<ref>)
for ψ(x,0).

We start with the solution for x>d (<ref>).
If we rewrite K(x,p,t)  for q^2=p^2-2 m V < 0 and use 

    X_l≡ d(2l+1) ,

 we find  (see section 4.4 in <cit.> for details)


    K(X_l,p,t)=U_0+U_1-U_2 ,    
        U_0=e^-i V t/ħ1/√(2 πħ)
     e^-i t q^2/2 ħ m+i X_l q/ħ =
     1/√(2 πħ)
     e^- i p^2 t/2 m ħ-X_l√(2 m V-p^2)/ħ ,
    
        U_1=
     κ/√(2 π t  i)1/√(2 πħ)
    e^-i V t/ħ∫_0^∞  e^i (X_l+u)^2 κ ^2/2 t
    e^i √(2 m V-p^2) u/ħ
    du , 
       U_2=κ/√(2 π t  i)1/√(2 πħ)
    e^-i V t/ħ∫_0^∞  e^i (X_l-u)^2 κ ^2/2 t
    e^i √(2 m V-p^2) u/ħdu
      .

Inserting 

    κ/√(2 π t  i)e^i (X_l± u)^2 κ ^2/2 t=
      ∫ _-∞^∞1/2 πħ
      e^-i q^2 t/2 m ħ+i q X_l/ħ±i q u/ħdq ,

we obtain

    U_1-U_2=
     1/(2 πħ)^3/2e^-i V t/ħ∫_-∞^∞ 
     e^-i q^2 t/2 m ħ+i q X_l/ħ2 i ħ q/2 m V-p^2+q^2dq  .


 We find for  U_1-U_2, that contains no oscillatory part
 
    ∫_-∞^∞(U_1-U_2)f(x-d,p)dp≈ 0 ,

 since
 
    ∫_-∞^∞f(p) p^n e^i p (x-d)/ħ dp= 
      (-i ħ∂/∂ x )^nψ(x-d,0)=0      x > d , n=0,1,2,N ,

 and so the contribution to (<ref>) will be proportional to
 Δ p^N+1 ,
 if f(p) decays rapidly enough so that the integrals
 
    ∫_-∞^∞f(p) p^n dp

  exist up to n=N+1.
 

 If ψ(x) is infinitely differentiable for
 x>d, th eleft hand side of (<ref>) vanishes identically.
So we conclude that for wavepackets that are sufficiently smooth for x > d , U_0 is the only relevant contribution.

We obtain for the wavefunction for x > d


    ψ(x,t) ≈∑_l=0^∞∫_0^t1/√(2 πħ)
     e^-X_l√(2 m V-p_0^2)/ħ∫_-∞^∞e^- i p^2 (t-τ)/2 m ħ
     e^i p (x-d)/ħ
     f(p) dp · g_l(τ)dτ .


This result contains the free time evolution  of the initial wave packet
that is shifted by d to the right. Proceeding as for the asymmetric square well(see <cit.>) we can approximate the 
convolution integral by

    ∫_0^t1/√(2 πħ)∫_-∞^∞e^- i p^2 (t-τ)/2 m ħ
     e^i p (x-d)/ħ
     f(p) dp · g_l(τ)dτ≈
       ∫_0^∞1/√(2 πħ)∫_-∞^∞e^- i p^2 (t-τ)/2 m ħ
     e^i p (x-d)/ħ
     f(p) dp · g_l(τ)dτ≈
       1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^i p (x-d)/ħ R(p)^2l (1-R(p)^2)f(p) dp 
      
       ≈ R(p_0)^2l (1-R(p_0)^2) ·1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^i p (x-d)/ħf(p) dp  ,

where R(p) is given by 

    R(p)=-1 + 2 k-2 √(k(k-1))     k=p^2/2 m V .

Here the assumptions about the concentration
of the wave packet (<ref>, <ref>) justifies putting R(p) before 
the integral. Evaluating the sum in (<ref>), we obtain

    ψ(x,t) ≈∑_l=0^∞∫_0^t1/√(2 πħ)
     e^-X_l√(2 m V-p_0^2)/ħ∫_-∞^∞e^- i p^2 (t-τ)/2 m ħ
     e^i p (x-d)/ħ
     f(p) dp · g_l(τ)dτ≈
    e^- d √(2 m V-p_0^2)/ħ/1-e^-2 d √(2 m V-p_0^2)/ħR^2(p_0)
     (1-R^2(p_0))
     1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^i p (x-d)/ħf(p) dp  .

For the neglected part of the time integral (<ref>)

    Δ_l≡∫_t^∞1/√(2 πħ)∫_-∞^∞e^- i p^2 (t-τ)/2 m ħ
     e^i p (x-d)/ħ
     f(p) dp · g_l(τ)dτ ,

we find the estimates (see Appendix <ref>)

    |Δ_l|∼ O(l)· 
    O([2 ħ/V t]^1/4)

Since 

    ∑_l=0^∞ 
    e^-X_l√(2 m V-p_0^2)/ħΔ_l∼e^d √(2 m V-p_0^2)/ħ/(-1+e^2 d √(2 m V-p_0^2)/ħ)^2· 
    O([2 ħ/V t]^1/4)

the approximation for ψ(x,t) given by (<ref>)
can be used for t ≫ħ/V. 

The wavepacket <ref> describes a free evolving wavepacket that started at t=0 at the position x_0+ d. It leaves the 
barrier without any time delay and is instantaneously transmitted through the barrier attenuated 
by a factor of the magnitude

    e^ -d √(2 m V-p_0^2)/ħ .



So within our approximations we find that the tunneling time is zero, where we assume in accordance with <cit.> 
that the transmission begins when a freely evolving 
wave packet starting at t=0 at the position x_0 would have arrived at the barrier, or equivalently when a classical
particle starting at t=0, x_0=0 has reached the barrier. It is not practicable to follow the peak of the original wavepacket until the 
entrance of the tunnel since this peak will be deformed by oscillations during the reflection process(<cit.>,<cit.>) and also the position expectation value
might undergo a slight attenuation before the tunnel as it was shown for infinite walls (<cit.>).


Performing the inverse Laplace transform for x<0 and 0 < x < d 
(<ref>,<ref>)and making the same approximations 
as in the previous case we find for the wavepacket
to the left of and within the barrier:


    ψ (x,t)≈κ/√(2 π t  i)∫_-∞^0  e^i (x-y)^2 κ ^2/2 tψ(y,0) dy
    +
      1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^-i p x/ħ R(p)f(p) dp + 
       ∑_l=0^∞
     e^-2 d (l+1) √(2 m V-p_0^2)/ħ1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^-i p x/ħ R(p)^2 l +1(R(p)^2-1)f(p) dp ≈
       κ/√(2 π t  i)∫_-∞^0  e^i (x-y)^2 κ ^2/2 tψ(y,0) dy +  
       {1+e^- 2 d √(2 m V-p_0^2)/ħ/1-e^-2 d √(2 m V-p_0^2)/ħR(p_0)^2·
    (R^2(p_0)-1)
     }
      R(p_0)1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     e^-i p x/ħ f(p) dp 
        x < 0  ,


    ψ (x,t)≈∑_l=0^∞
     e^-(2 d l+x) √(2 m V-p_0^2)/ħ1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
      R(p)^2 l(R(p)+1)f(p) dp ≈
        -∑_l=0^∞
     e^-(2 d (l+1)+x) √(2 m V-p_0^2)/ħ1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
     R(p)^2 l+1(R(p)+1)f(p) dp ≈
       {1/1-e^-2 d √(2 m V-p_0^2)/ħR(p_0)^2 
    - e^- 2( d-x) √(2 m V-p_0^2)/ħ/1-e^-2 d √(2 m V-p_0^2)/ħR(p_0)^2·
     R(p_0)
     }
        · (R(p_0)+1) e^- x√(2 m V-p_0^2)/ħ· 1/√(2 πħ)∫_-∞^∞e^- i p^2 t/2 m ħ
      f(p) dp 
        0 < x < d  .


The solution within the barrier differs from the solution of the potential step (see <cit.>) only
by a time-independent factor


    ψ_step(x,t)≈ (1+R(p_0))
    e^-√(2 m V-p_0^2) x/ħ1/√(2 πħ)∫_-∞^∞  e^- i p^2 t/2 m ħ
    f(p) dp .


So we see that the time it
needs until the barrier is (approximately) empty again, is
independent of the thickness of the barrier.
The solution on the left side consists of the incoming and the reflected wave packet (<ref>).
In contrast to the potential step the reflected wavepacket experiences a permanent
attenuation, since a part of the wavefunction has tunneled through the barrier.
After the reflection process the wavefunction consists of a reflected and a transmitted
wavepacket (<ref>) only.
An explicit calculation yields


    | {1+e^- 2 d √(2 m V-p_0^2)/ħ/1-e^-2 d √(2 m V-p_0^2)/ħR(p_0)^2·
    (R^2(p_0)-1)
     }
      R(p_0) |^2+ 
       |e^- d √(2 m V-p_0^2)/ħ/1-e^-2 d √(2 m V-p_0^2)/ħR^2(p_0)
     (1-R^2(p_0))|^2 = 1  ,


which confirms that the integral over the probability density is conserved and our approximations
are consistent.


§ THE TUNNELING TIME


Within our approximations we found the tunneling time to be zero, 
since our solution (<ref>)
indicates that the wavepacket leaves the tunnel, shifted by d to the right. 
Here we have assumed, 
that all functions of R(p) can be pulled out of the integral (<ref>), and 
therefore 
they do not influence the 
dynamics of the wave packet. If we take into account the first order contributions of R(p), 
we find
a small, but finite tunneling time. 
We will from now on assume that the initial wave function
is an uncorrelated function of the form

    f(p)=G(p)e^-i p x_0/ħ ,

where G(p) is a real function that yields a momentum expectation value p_0. The position expectation value is then given by x_0.

The impact of an additional factor Z(p) on the initial wavepacket f(p) is twofold

    Z(p)f(p)=z(p)e^i μ (p) f(p) ,   z(p)=|Z(p)| .


The probability density in momentum space is only affected by the absolute value z(p).
If we restrict ourselves to first order contributions in p-p_0, 

    z(p) ≈ z(p_0)+ z_1(p-p_0)


we find that we can neglect
the momentum shift:

    N^2≡∫_-∞^∞|f(p)|^2 z(p)^2 dp
     ≈  z(p_0)^2∫_-∞^∞|f(p)|^2 dp =
     z(p_0)^2 
       ⟨p̂⟩ = 
    N^-2∫_-∞^∞p|f(p)|^2 z(p)^2 dp 
       ≈
    p_0+2 (z(p_0 ))^-1 z_1∫_-∞^∞|f(p)|^2 p(p-p_0) dp ≈
    p_0 .

 Using also a linear approximation  for μ'(p)
 
    μ'(p) ≈μ_1+μ_2(p-p_0)

we find for the position expectation value
 
    x_0 =
      N^-2∫_-∞^∞ z(p) f^*(p)e^-i μ (p)(i ħ∂/∂ p)
      z(p) f(p)e^i μ (p) dp = 
       N^-2 x_0∫_-∞^∞G(p)^2 z(p)^2 dp
     -N^-2ħ∫_-∞^∞G(p)^2 z(p)^2 μ'(p)dp
    ≈    x_0-ħμ_1 .


Therefore the phase of the functions

    R(p)^2l (1-R(p)^2)

in (<ref>)
yields a shift of the position of each particular wave packet constituting  (<ref>)

    ψ(x,t) ≈∑_l=0^∞1/√(2 πħ)e^-X_l√(2 m V-p_0^2)/ħ∫_-∞^∞e^- i p^2 t/2 m ħ
     e^i p (x-d)/ħ R(p)^2l (1-R(p)^2)f(p) dp  .


We find

    Arg[1-R(p)^2]=
    Arctan[-1+2 l/2 √(k-k^2)] ≈ 
    Arctan[-1+2 k_0/2 √(k_0-k_0^2)]+
    2/√(2 m V-p^2_0)·(p-p_0) 
       Arg[R(p)]=
    -Arccos[1-2 k] ≈ 
    -Arccos[1-2 k_0]+
    2/√(2 m V-p^2_0)·(p-p_0)  ,

where we have used the definitions

    k=p^2/2 m V ,   k_0=p_0^2/2 m V .

So we see that the lth term of (<ref>) will experience a phase shift of 

    2 (1+2l) (p-p_0)/√(2 m V-p^2_0) ,

corresponding to a translation by

    δ x_l = 2 (1+2l) ħ/√(2 m V-p^2_0)=2 (1+2l) ħ/p_0√(k_0/1-k_0)


to the left. So the delay time will be

    T_l=  2 (1+2l) m ħ/p_0√(2 m V-p^2_0)

instead of zero.

Each term is attenuated by a factor of the magnitude

    e^ -X_l√(2 m V-p_0^2)/ħ .

For thick barriers, if

    d √(2 m V-p_0^2)/ħ≫ 1 ,


the first term will dominate the sum (<ref>), and the tunneling time will be given by

    T_0=2 m ħ/√(2 m V-p^2_0)p_0


This is exact the time Hartmann <cit.>found as upper limit for the tunneling 
time through
thick barriers.
If the barrier gets thinner the other wavepackets for l>0
will also come into play (see figure <ref>). Each of them leaves the barrier at a different time T_l.
But since they appear very shortly after each other they may appear as one smeared out wavepacket
with a delay time bigger than T_0.

 Note that apart from the absolute and relative magnitude of the wavepackets, all characteristic
quantities of tunneling are independent of the width d. In the case of the wavefunction pictured in figure <ref> where the uncertainties are
given by 

    Δ x≈√(a/2) ,  Δ p ≈√(ħ/2 a) ,

the ratio between the distance of the centre of the wavepackets and the position uncertainty the determines the distinguishability between the wavepackets reads

    δ x_l+1-δ x_l/Δ x= 
     4 ħ/p_0Δ x√(k_0/1-k_0)≈8 Δ p/p_0√(k_0/1-k_0) .









§ DISCUSSION AND CONCLUSIONS 

Our result for the  tunneling time (<ref>) is not an exact reproduction
of Hartmann's result (<cit.>,<cit.>) which predicts an increasing tunneling time with the thickness of the barrier 
before saturation 
takes place. However Spielmann et al. found a decreasing tunneling time in their experiments with
electromagnetic waves propagating through photonic band gap materials (see <cit.> fig. 3, <cit.> fig 1 ).
This qualitative behaviour is in good agreement with our results that also predict a decreasing tunneling time since for
thicker barriers the later wave packets are more and more attenuated.







For thicker barriers the conclusion of both calculations is that the tunneling time for sufficiently
peaked wave packets is given by (<ref>).
This is as far interesting as the results were obtained by completely different methods. Moreover we ensured in our calculations
that the initial wavepacket is only located at the left side of the barrier (<ref>) which is not clearly
guaranteed by Hartmann's approach. So our result makes sure that the Hartmann time (<ref>) is not some relic of the parts of 
the initial wavepacket that were at the right hand side
of the barrier from the beginning. 

We did not take into account the parts of the initial wavepacket with energies 
near or greater than the critical energy of standard transmission through the barrier
(see <ref>). In Appendix <ref> we have derived a consistency condition for this 
approximation for the case of a special class of initial wavepackets. So this is a further counterexample to the idea that 
only energies greater than the barrier height contribute to tunneling (<cit.>).
                      


We also found out that the approximate solution within a finite barrier differs
from the solution within the potential step only by a time-independent factor (<ref>,<ref>)
which also indicates that important
dynamical properties are independent of the thickness of the barrier.
It would be especially interesting if
this is also true for more general tunneling processes as the tunneling out of a potential well that could model
radioactive decay or tunneling out of atoms as provided by the attoclock experiment  <cit.>. Moreover an application of
the method of Laplace transformation to relativistic wave equations would yield a picture of the reflection and tunneling processes
in the relativistic case.


§ ACKNOWLEDGMENTS


We thank the referees for their comments and suggestions that helped to improve this article.


§ ESTIMATION OF THE CONVOLUTION INTEGRAL
 


According to <cit.>, the inverse Laplace transform of ρ(s)^l is given by 

    ℒ^-1( ρ(s)^l)= 
    l/i^l tJ_l[V t/2 ħ]e^-i V t/2 ħ .


Therefore we find

    g_l(t)=ℒ^-1(ρ^(2 l)(s)(1-ρ^2(s)))=
    l/i^(2l) tJ_2l[V t/2 ħ]e^-i V t/2 ħ-
    l/i^(2(l+1)) tJ_2(l+1)[V t/2 ħ]e^-i V t/2 ħ .


For an integral of the form

    u(l,t)=
    ∫_t^∞ e^∓ i p Q/ħ
    e^-i p^2 (t-τ)/2m ħ e^-i V τ/2 ħl/i^l  τ
    J_l[V τ/2 ħ] dτ ,



we get the following estimate


    |u(l,t)|=|
    ∫_t^∞e^i p^2 τ/2m ħ-i V τ/2 ħl/τ
    J_l[V t/2 ħ] dτ| ≤∫_t V/2 ħ^∞|l J_l(y)/y| dy ≤
       (∫_t V/2 ħ^∞l/y^1+2ϵ dy )^1/2·(∫_0^∞(J_l(y))^2/y^1-2ϵ)^1/2 = 
        l 1/√(2 ϵ)(2 ħ/V t)^ϵ·(
    2^2 ϵΓ [1-2 ϵ]Γ[l+ϵ]/2 (Γ[1-ϵ])^2 Γ[1+l-ϵ])^1/2 ≤
        l 1/√(2 ϵ)(2 ħ/V t)^ϵ·(
    2^2 ϵΓ [1-2 ϵ]/2 (Γ[1-ϵ])^2 )^1/2     0 < ϵ < 1/2 ,



where we have applied the Schwarz inequality and the integral formula
(<cit.>)

    ∫_0^∞(J_l(y))^2/y^1-2ϵ=
    2^2 ϵΓ [1-2 ϵ]Γ[l+ϵ]/2 (Γ[1-ϵ])^2 Γ[1+l-ϵ] .


So we conclude setting ϵ=1/4 for Δ_l (<ref>)

    |Δ_l|∼ O(l)· 
    O([2 ħ/V t]^1/4)




§ DECREASING BEHAVIOUR IN MOMENTUM SPACE OF FUNCTIONS WITH COMPACT SUPPORT AND THE EXAMINATION OF THE USED
ASSUMPTIONS ABOUT ORDERS OF MAGNITUDE BASED ON AN EXAMPLE
 

According to the Palay Wiener theorem (<cit.>), functions with compact support in position space can not be restricted 
to a finite interval in momentum space as well.
 Nevertheless the concentration of those functions in momentum space around their 
expectation value can be shown 
to be prescribed by the Fourier transform of a  generic  reference function. We choose an appropriate wavepacket with a reference 
function of Gaussian shape and derive the conditions that justify (<ref>,<ref>,<ref>). We also show that we can use the reference function
for the evaluation of (<ref>) as we did for the example presented in figure <ref>.


We start with a normalized wavefuntion ψ(x) with position expectation value x_0 that 
is assumed to be zero for |x-x_0| ≥ B. Let ψ_0(x) be a generic reference function that fulfills

    [    ψ_0(x)=ψ(x)                   |x-x_0| ≤ B;  ψ_0(x)=-δψ(x)                |x-x_0| > B  , ]
  
where δψ is zero for |x-x_0| ≤ B. Then ψ(x) is represented by the sum

    ψ(x)=ψ_0(x)+δψ (x)



If the norm of δψ is given by ϵ we find for the reference function


    ⟨ψ_0 |ψ_0⟩ =1+ϵ^2 .



The representation of ψ(x) in momentum space reads


    f(p)=f_0(p)+δ f(p)

where f_0(p) and δ f(p) are the representations  ψ_0(x) and δψ(x) in momentum space.

Using the Schwarz inequality and taking into account that the Fourier transform preserves the L_2 norm we 
get for the integral over the momentum density in the region outside the interval (p_min,p_max)


    ∫_-∞^p_min f^*(p)f(p)  dp + ∫_p_max^∞ f^*(p)f(p)  dp  ≤
       ∫_-∞^p_min f_0^*(p)f_0(p) dp+ ∫_p_max^∞ f_0^*(p)f_0(p) dp +ϵ^2 +
    2 √(1+ϵ^2)ϵ .




In order to provide an explicit example, we choose the initial wavepacket with the position and momentum
expectation values x_0,p_0



    ψ (x)=N^-1/2 
    e^-i p_0 x/ħ
     e^-(x - x_0)^2/ 2 a[-L + x - x_0/√(a)]^2 
    [L + x - x_0/√(a)]^2          |x-x_0| < L √(a)
       ψ(x)=0         |x-x_0| ≥ L √(a) ,

 
 where the normalization constant N  is given by
 
    N= √(a)/16 e^-L^2 (2 L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
       e^L^2 (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L]) .

The factors multiplied to the Gaussian ensure that ψ(x) is continuously differentiable at x=x_0± L √(a).

The position and momentum uncertainty are determined by



    (Δ x) ^2  =
    a (2 L (-945 + 210 L^2 - 52 L^4 + 8 L^6) + 
     e^L^2 (945 - 840 L^2 + 360 L^4 - 96 L^6 + 16 L^8) √(π)[L])/
    4 L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
     2 e^L^2 (105 - 120 L^2 + 72 L^4 - 32 L^6 + 16 L^8) √(π)[L]
       (Δ p) ^2 =
    ħ^2 2 L (-225 + 18 L^2 + 12 L^4 + 8 L^6) + 
     e^L^2 (225 + 8 L^2 (-21 + 5 L^2 + 4 L^4 + 2 L^6)) √(π)[L])/2 a (2 L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
       e^L^2 (105 - 120 L^2 + 72 L^4 - 32 L^6 + 16 L^8) √(π)[L]) ,



where

    lim_L →∞(Δ x) ^2 =a/2       Δ x ^2≤a/2     L > 2
      
       lim_L →∞(Δ p) ^2 =ħ^2 /2 a       Δ p ^2≤ 1.2 ħ^2/a     L > 2


So we see that it is possible to choose the momentum uncertainty sufficiently small to justify (<ref>).
We take as  reference function ψ_0(x) the extension of (<ref>) to the whole real line. We find for the norm of the difference
function δψ according to (<ref>)

    ϵ^2 =
    2 L (105 - 50 L^2 + 20 L^4 - 8 L^6) + 
     e^L^2 (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L]
      /
    2 L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
     e^L^2 (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L] ,


where                                                        
 
    lim_L →∞ϵ =√(24/√(π))e^-L^2/2L^-9/2  ϵ≤√(24/√(π))e^-L^2/2L^-9/2     L > 3

For the further calculations we introduce the dimensionless quantities

    P=√(a)/ħ p ,     P_0=√(a)/ħ p_0 ,
        X=x/√(a) ,    X_0=x_0/√(a) ,   D=d/√(a) .

We find for the momentum representation of the reference function ψ_0(x) 

    f(p)=a^1/4/√(h)4  e^-(1/2) (P - P_0) (P - P_0 + 2 I X_0) (3 + L^4 + 
       2 L^2 (-1 + (P - P_0)^2) - 6 (P - P_0)^2 + (P - P_0)^4)/√(
    2 e^-L^2  L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
      (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L]) .

Assuming that the interval (p_min,p_max) is symmetric around p_0 > 0 we find
for the probability outside the interval

    𝒫_rest≡∫_-∞^p_min f_0^*(p)f_0(p)+
     ∫_p_max^∞ f_0^*(p)f_0(p) dp =
     
       {
     2 e^-P_0^2(1+(K-2)K) 
    (-1 + K) P_0 (-39 + 72 L^2 - 88 L^4 + 32 L^6 +. .
       . 2 (-1 + K)^2 (83 + 24 L^2 (-3 + L^2)) P_0^2 + 
        4 (-1 + K)^4 (-17 + 8 L^2) P_0^4 + 8 (-1 + K)^6 P_0^6 ) 
     
        
      +.
     (105 + 
          8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[P_0 (K-1)]
      }·
       (
     Abs[
    2 e^-L^2  L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
       (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L]
     ]
        )^-1 ,

where K ≡ p_max/p_0.

According to (<ref>) the application of the approximation (<ref>) means that we neglect a portion of the probability 
of the magnitude ϵ^2 + 𝒫_Rest. Since the wavepackets that leave the barrier are of the magnitude 
 
    e^ -d(2l+1) √(2 m V-p_0^2)/ħ=e^-(2l+1) D P_0√(1-k_0/k_0) .

  our results are relevant compared to the neglected parts if the condition

  
    D P_0 (2l+1)
     √(1-k_0/k_0)≲[ϵ+𝒫_Rest]

  
  is fulfilled.

For the initial wavfunction evaluated in fig. <ref> with the parameters 

    P_0=10, k_0=1/2, X_0=-20, L=20

and with the choice 

    K=p_max/p_0=1.4 ,

(<ref>) reads

    D (2l+1)· 10  ≲  18 .


So the wavepackets up to  D (2l+1)≤ 1.8 meet <ref>. Moreover the set of parameters fulfills the requirements of (<ref>), since

    1/√(1-p_max^2/2 m V)=
    1/√(1-p_max^2/p_0^2)≈ 7 .


For the evaluation of the first three terms of (<ref>) in fig.<ref> we have used f_0(p) instead of f(p): This is justified since for our choice
of parameters (<ref>)  ϵ≪ 1 . Moreover the obtained position shift (<ref>)

    δ x_l = 2 (1+2l) ħ/√(2 m V-p^2_0)= 
    2 (1+2l)√(a)/P_0√(k_0/1-k_0)

is bigger than the correction of the positions expectation value x_0 caused by δψ. We find


    δ x_ref≡⟨ψ_0  x ψ_0⟩ - ⟨ψ  x ψ⟩ = ⟨δψ  x δψ⟩ =
     
        x_0·2 L (105 - 50 L^2 + 20 L^4 - 8 L^6) + 
     e^L^2 (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L]/
    2 L (-105 + 50 L^2 - 20 L^4 + 8 L^6) + 
     e^L^2 (105 + 8 L^2 (-15 + 9 L^2 - 4 L^4 + 2 L^6)) √(π)[L] .

Since 

    | δ x_ref| ≤|x_0|·24/√(π)e^-L^2L^-9     L > 2 ,

the condition δ x_ref ≪ δ_l is fulfilled if

    24/√(π)e^-L^2L^-9 ≪ 2(1+2l) P_0^-1|X_0^-1| 
     √(k_0/1-k_0) ,


which is the case for the choice of parameters (<ref>).





99
HaugeE.H.Hauge,J.A. Stovneng (1989):
Tunneling Times: A critical review, Rev.Mod.Phys. 61, 917
Landauer  R.Landauer , Th.Martin (1994): Barrier interaction time in tunneling, Rev.Mod.Phys.66, 217
Winful H.Winful (2006): Tunneling time, the Hartmann effect and superluminality:
 A proposed resolution of an old paradox, Phys. Rep.436,1-69
 Olkhovsky V. Olkhovsky et. al. (2004): Unified analysis of photon
 and particle tunneling Phys. Rep. 398, 133-178
 Razavy M.Razavy: Quantum theory of Tunneling, World Scientific Singapore 2003
  GalaponE.Galapon (2012): Only above barrier energy components contribute
  to barrier traversal time, Phys.Rev.Lett.108,170402
 Landsmann A. Landsman, U.Keller (2014: Tunneling time in strong field ionization, J.Phys.B 47, 204024
 Landsmann1 A. Landsman, U.Keller (2015: Attosecond science and the tunneling time problem, Phys.Rep.547,1-24
Wang Z.S.Wang et.al.(2004: Quantum tunneling time, Phys.Rev. A 69, 052108
Cheng Y. Cheng (2010): On the tunneling time of arbitrary continuous potentials and the Hartmann effect, Chin.Phys.B, 19,No.11, 117305
Pfeiffer A.Pfeiffer et. at. (2012):
Attoclock reveals natural coordinates of the 
laser-induced tunneling current flow in atoms, Nature Physics Vol. 8, 76
Eckle P.Eckle et.al. (2008): Science Vol.322 1525
Landsmann2 A.Landsman (2014): Ultrafast resolution
Hartmann T. Hartmann(1962): Tunneling of a Wave Packet, J. Appl. Phys. 33, 3427
of tunneling delay time, Optica Vol.1,No.3,343
  MacColl L.A.MacColl(1932): Note on the Transmission and Reflection of Wave Packets by Potential
                  Barriers, Phys.Rev. 40/ May 1932/621-626
 
 Hasan M. Hasan, B.Mandal (2018): Phys.Lett. A382, 248-252
 SteinbergA. Steinberg et. al. (1993): Measurement of single-photon tunneling time, Phys.Rev.Lett. 71/5,708
 Spielmann Ch. Spielmann, R. Szipös, A.Stingl, F. Krausz (1994) :
Tunneling of optical pulses throug photonic band gaps, Phys.Rev. Lett. 73, 2308
Longhi S.Longhi et. al (2001):Superluminal optical pulse propagation, Phys. Rev. E 64055602
Chiao R.Y. Chiao et. al.(1993): Faster than light?, Scientific American 269, 52
Winful1 H.G. Winful (2003): Mechanism for 'superluminal' tunneling, Nature 424,638

 Methode N.Riahi (2017): Solving the time-dependent Schrödinger
equation via Laplace transform, Quantum Stud.: Math. Found.  4/2, 103-126

Bernardini A.E. Bernardini (2009): Stationary phase method and delay times for relativistic
 and non-relativistic tunneling particles, Ann.of Phys. 324, 1303-1339
Pereyra P.Pereyra (2000): Closed formulas for tunneling time in superlattices,
Phys. Rev.Lett. 84/8, 1772
LosTunnel
V.F.Los,M.V. Los (2012) : A time-dependent exact solution for wave-packet scattering at a rectangular barrier,
J.Phys.A: Math.Theor. 45,095302
 Bounce
 M.A.Doncheski et. al. (1999): Anatomy of a quantum 'bounce',Eur.J.Phys.20, 29-37
 Erdely A.Erdely: Tables of Integral Transforms, McGraw-Hill Book Company, 1954
 Magnus W.Magnus, F.Oberhettinger, R.P.Soni: Formulas and Theorems for the 
 Special Functions of Mathematical Physics, Springer-Verlag 
  Za A.Zayed: Handbook of Function and Generalized Function Transformations, CRC Press 1996 


 
 


