V Gradinaru

G A Hagedorn

A Joye

Tunneling dynamics and spawning with adaptive semi-classical wave-packets Tunneling dynamics and spawning with adaptive semi-classical wave-packets

Keywords: numbers: 03.65.Sq, 82.20.Wt, 02.70.Hm, 02.60.Cb semi-classical, time-dependent Schrödinger equation, wave-packets, tunneling, spawning, timedependent discrete variable representation

Tunneling through a one-dimensional Eckart barrier is investigated using a recently developed propagation scheme based on semi-classical wave-packets. This version of the time-dependent discrete variable representation method yields linear equations for the parameters, is fully adaptive, and does not require a frozen Ansatz in order to approximate the exact solution of the Schrödinger equation accurately. We rely on an analytical result to derive a new algorithm to spawn a second family of semi-classical wave-packets after the tunneling has occurred. Numerical results for a benchmark problem demonstrate the accuracy of the new method.

I. INTRODUCTION

Semi-classical wave-packets were introduced [START_REF] Hagedorn | [END_REF]2 in order to deal with the time-dependent Schrödinger equation in the semi-classical scaling, i.e.,

i ε ∂ψ ∂t = Hψ, (1)
where ψ = ψ(x, t) is the wave-function depending on the spatial variables x = (x 1 , . . . , x N ) ∈ R N and the time t ∈ R. The Hamiltonian operator H, which depends on ε, is H = T + V with the kinetic and potential energy operators

T = - N j=1 ε 2 ∂ 2 ∂x 2 j and V = V (x),
where the real-valued potential V acts as a multiplication operator on ψ in L 2 (R N ).

In molecular quantum dynamics, (1) is a Schrödinger equation for the nuclei on an electronic energy surface in the time-dependent Born-Oppenheimer approximation [START_REF] Hagedorn | Spectral theory and mathematical physics: a Festschrift in honor of Barry Simon's 60th birthday[END_REF][4][5][START_REF] Tannor | Introduction to Quantum Mechanics: A Time Dependent Perspective[END_REF][START_REF] Teufel | Adiabatic perturbation theory in quantum dynamics[END_REF] . In this situation, ε 2 is the mass ratio 8 between the electrons and nuclei, of magnitude 10 -4 .

Semi-classical wave-packets were employed for an analytic proof that the wave function can be approximated with high asymptotic accuracy in ε by complex Gaussians times polynomials [START_REF] Hagedorn | [END_REF]2 . Recently, they have been put to use in a numerical scheme for multi-particle quantum dynamics in the semi-classical regime 9 . In the quadrature version of this algorithm, the relation with the timedependent discrete variable representation TDDVR 10,11 and quantum-dressed classical mechanics 12 is evident. The wave-packets constitute an orthogonal basis that adapts in time to the evolution at hand. Time-dependent adaptive spaces and grids arise in the same natural manner. In one space dimension, the semi-classical wave-packets are just scaled and shifted Hermite polynomials times complex Gaussians, but in higher space dimensions they are both more general and more suitable than tensor products of Hermite functions 9 . An advantage of the semi-classical wave-packets is the freedom given by their supplementary parameters. This freedom yields linear equations for parameters that are fully adaptive and do not have to be fixed as sometimes happens in the TDDVR 10 . Making use of this freedom, the algorithm in Ref. 9 proves not only suitable for higher-dimensional cases, but it is time-reversible and ensures an unitary propagation, hence perfect norm conservation. As long as the approximation remains valid we even have no drift in the energy. Last but not least: the classical picture of the dynamics is ε-blurred: for ε → 0 we get classical dynamics with its advantageous numerical propagator: Störmer-Verlet. Larger ǫ puts in more and more quantum effects, and it complies well with the Born-Oppenheimer-Approximation.

Methods based on wave-packets suffer if the widths of the basis functions become too large. Tunneling thus seems difficult to address, since the occurrence of a delocalized wave function would imply large width and require a large number of wave-packets 10,11 . Tunneling through the Eckart potential is, on the one hand, a simple model for bi-molecular reaction dynamics (e.g. H + H 2 exchange reaction) and on the other hand, a non-trivial benchmark test for the accuracy of the semi-classical approximation 10,11,13,14 . So, we stick to the one dimensional case in this work and use the Eckart potential for the numerical experiments.

Let us anticipate here the results in the notation of the next sections that explain the technical details. Figures 1 and2 show the squared absolute value of the wavefunction together with its representation in terms of coefficients c k of basis functions at the some times shortly after the tunneling. The initial state for Figure 1 was a typical Gaussian. Note the difference between scales of the reflected and tunneled parts of the wave-function. The higher order coefficients are responsible for the tunneled part, but not all 512 considered terms are necessary. The dashed line that appears on the right domain represents the spawned Gaussian as computed by our algorithm, whereas the dotted line displays the potential. Figure 2 reflects the case when we take as an initial wave- packet a polynomial times a Gaussian with same position and width parameters as in the previous case but slightly less momentum, in order to have only slightly less total energy in the system. The results show that our spawning technique reduces dramatically the number of used basis functions, while keeping the error resonably small. More numerical results and the exact values of the parameters used are presented in the last section.

The main difficulty in the simulation is the fact that typical solutions have wavelength of order ǫ, move with velocity of order 1 in a non-confining potential, and have spreading width starting from order √ ǫ. Note that too small ǫ and large final times make the Fourier method computationaly intensive. The existent wave-packets methods behave better for smaller ǫ but need many basisfunctions for a spreading wavefunction as in the tunneling problem. In analogy with the multiple spawning approach [START_REF] Ben-Nun | [END_REF] , we give a rigorously justified spawning algorithm which specifies when and how to add a new family of semi-classical wave-packets during the tunneling in one dimension.

In the next section we briefly discuss the semi-classical wave-packets and sketch their use for a practical solution of the problem (1). We make use of a mathematically rigorous proof 17 that in the 1D tunneling problem with a smooth bounded potential, the part of the wave function that tunnels is a Gaussian to leading order in ε at sufficiently large times. This analytical result gives formulas for the emerging wave function that unfortunately are not of practical use. The question of how to spawn the new Gaussian is answered in the third section by an algorithm involving the semi-classical wave-packets. The same algorithm is valid for the case of a reflected Gaussian. Numerical experiments validating our new method are presented in the last section.

II. THE ANSATZ A. Semi-classical wave-packets

Since the analytical tunneling result has been proved only in the one dimensional case, we stick here to N = 1 and refer to Ref. 2 and 9 for the construction and propagation algorithms in the higher dimensional case. A Gaussian wave-packet is parametrized as 18

ϕ ε 0 [q, p, Q, P ](x) = (πε) -1/4 (Q) -1/2 • exp i 2ε P Q -1 (x -q) 2 + i ε p(x -q) , (2)
where q, p ∈ R represent the position and momentum, respectively, and Q, P ∈ C satisfy the compatibility conditions:

QP -P Q = 0, (3) QP -P Q = 2i . (4)
The last two relations are equivalent to requiring that

Y = Re Q Im Q Re P Im P
be symplectic, an important property for the numerical integration 19 of conservation laws. A complete L 2 -orthonormal set of functions

ϕ k (x) = ϕ ε k [q, p, Q, P ](x) (5)
for non-negative integers k, is recursively constructed as follows 2 : Let x denote the position operator (acting on functions of x by multiplication with x) and y = -iε∇ x the momentum operator, and introduce the raising operator R and lowering operator L as

R = - i √ 2ε P (x -q) + Q(y -p) , L = i √ 2ε (P (x -q) + Q(y -p)) .
Define

ϕ k+1 = 1 √ k + 1 R ϕ k . (6)
It then turns out that these functions are orthonormal, as the eigenfunctions of the hermitian operator LR = RL + I. Moreover, we have

ϕ k-1 = 1 √ k L ϕ k ,
(the right-hand side is zero if k = 0), and the functions ϕ k are polynomials of degree k multiplied by the Gaussian ϕ 0 . Since the above relations imply (see Ref. 2, equation (3.28))

x

-q = ε 2 QR + QL , (7)
we obtain the recurrence relation

Q √ k + 1 ϕ k+1 (x) = 2 ε (x -q)ϕ k (x) -Q √ k ϕ k-1 (x) ,
which permits us to compute the functions ϕ k at any given value x. In the one dimensional case we get (scaled) Hermite polynomials times ϕ ε 0 [q, p, Q, P ]. Let us emphasize the differences between this and TD-DVR methods [10][11][12] : the role of a complex width parameter A(t) that is decomposed as Im A(t) and Re A(t) (and sometimes fixed as in the frozen Gaussian approach) in the TDDVR is played here by two complex parameters Q(t) and P (t) that will adapt themselves to the dynamics and keep the matrix Y (t) symplectic. Each state ϕ ε k [q(t), p(t), Q(t), P (t)] is concentrated in position near q(t) and in momentum near p(t); the position and momentum uncertainties are (k + [START_REF] Hagedorn | [END_REF] 2 )ε |Q(t)| and (k + [START_REF] Hagedorn | [END_REF] 2 )ε |P (t)|. We approximate solutions to the equation (1) in the form of a finite linear combination of wave-packets, with a common highly oscillatory phase factor,

u(x, t) = e iS(t)/ε K k=0 c k (t) ϕ ε k [q(t), p(t), Q(t), P (t)](x) .
(8) This Ansatz is motivated by the remarkable fact that in the case of a quadratic (possibly time-dependent) potential V , the functions e iS(t)/ε ϕ ε k [q(t), p(t), Q(t), P (t)] are exact solutions to the Schrödinger equation 2 if the position and momentum parameters follow the classical equations of motion q = p, ṗ = -∇V (q) ; the linearized equations of motion

Q = P, Ṗ = -∇ 2 V (q) Q ;
and S(t) = t 0 ( 1 2 p(s) 2 -V (q(s)) ds is the classical action. On the other hand, for a non-quadratic potential, the wave function can be expanded in the orthogonal basis of wave-packets with time-dependent coefficients, with parameters determined by the equations of motion corresponding to a local quadratic approximation of the potential.

B. Propagation

We now give a full algorithmic description of the adaptive time-stepping. Assume that the step-size ∆t is given, and let the real scalars q n , p n , S n , the complex Q n , P n , and the complex coefficient vector

c n = (c n k ) k∈K , K = {0, • • • , K}, be such that u n = e iS n /ε K k=0 c n k ϕ ε k [q n , p n , Q n , P n ]
is an approximation to the solution of the Schrödinger equation (1) at time t n = n∆t. To compute the approximation u n+1 at time t n+1 , we proceed as follows:

1. Compute q n+1/2 , Q n+1/2 , and S n+1/2,-via

q n+1/2 = q n + ∆t 2 p n , Q n+1/2 = Q n + ∆t 2 P n , (9)
S n+1/2,-= S n + ∆t 4 (p n ) 2 .
2. Compute p n+1 , P n+1 , and S n+1/2,+ via

p n+1 = p n -∆t ∇V q n+1/2 , P n+1 = P n -∆t ∇ 2 V q n+1/2 Q n+1/2 , ( 10
)
S n+1/2,+ = S n+1/2,--∆t V q n+1/2 .

Update the coefficient vector

c n+1 = (c n+1 k ) k∈K as c n+1 = exp -∆t i ε F n+1/2 c n , (11)
e.g., by a few steps of the Arnoldi iteration 9 . Here, F n+1/2 = (f kℓ ) k,ℓ∈K is the Hermitian matrix with entries

f kℓ = ϕ n+1/2 k W n+1/2 ϕ n+1/2 ℓ , (12)
where ϕ n+1/2 k = ϕ ε k [q n+1/2 , p n+1 , Q n+1/2 , P n+1 ] are the basis functions and

W n+1/2 (x) = V (x) -U n+1/2 (x)
is the remainder in the local quadratic approximation to V , given at q = q n+1/2 by U n+1/2 (x) = V (q) + ∇V (q) (x -q) + 1 2 ∇ 2 V (q) (x -q) 2 . Note that f kℓ actually depends only on q n+1/2 and Q n+1/2 , but not on p n+1 and P n+1 , since the imaginary parts in the arguments of the Gaussian cancel out in (12). 4. Compute q n+1 , Q n+1 , and S n+1 via

q n+1 = q n+1/2 + ∆t 2 p n+1 , Q n+1 = Q n+1/2 + ∆t 2 P n+1 , (13)
S n+1 = S n+1/2,+ + ∆t 4 (p n+1 ) 2 .
The algorithm is of second order accuracy in the parameters q, p, Q, P, S, and c k and enjoys a number of attractive conservation and limit properties: time-reversibility, symplecticity and L 2 -norm conservation, and robustness in classical limit ǫ → 0. Regarding the position and momentum parameters q and p, the algorithm coincides with the Störmer-Verlet method applied to the corresponding classical equations of motion. In the limit of taking the full basis set ϕ k , with all k ∈ N, the variational approximation used in the remainder propagator becomes exact and the algorithm converges towards the Strang splitting (or symmetric Lie-Trotter splitting) exp(-

i ε ∆tH) ≈ exp(-i ε ∆t 2 T ) exp(-i ε ∆tV ) exp(-i ε ∆t 2 T
) of the Schrödinger equation. We refer to Ref. 9 for more details and the treatment of the higher dimensional case that clarifies the flexibility of this Ansatz.

C. Tunneling

It is well known that a one-dimensional, semi-classical wave packet with low energy scattered on a potential barrier gives rise to an exponentially small tunneling piece, in the semiclassical limit, which propagates beyond the potential barrier V . Less well known is the fact that this tunneling piece always has a Gaussian shape, irrespective of the details of the incoming wave choice (8), in the semiclassical limit. This statement is true as soon as the tunneling piece is away from the potential barrier, not necessarily in a scattering regime, and the characteristics of the Gaussian are determined by the classical dynamics. The fact that Gaussians always emerge in semi-classical tunneling, which provides the theoretical foundation of the spawning, is proven in 17 . The result, in a setup suitable for the present analysis, reads as follows:

Assume the potential x → V (x) is an analytic function on R, such that lim |x|→∞ |V (x)| = 0, which has a unique maximum V 0 > 0 at x = 0. Let u(x, 0) be a semiclassical wave packet of the form (8), centered around (q, p) in phase space with q < 0 and p > 0, such that p 2 /2 < V 0 . Then, for t > t 0 > 0, the solution of (1) decomposes as u(x, t) = u ref (x, t) + u trans (x, t), where, for ε << 1, the reflected part, u ref (x, t), is negligible for x > 0 whereas u trans (x, t) is negligible for x < 0. Moreover, u trans (x, t) = Cε 3/4 e -γ/ε u G (x, t) + o(ε 3/4 e -γ/ε ), for some γ > 0 and C ∈ C, where u G (x, t) is a semi-classical Gaussian wave packet concentrated around a classical trajectory (q(t), p(t)). The other characteristics of u G can be computed and depend on the details, just as C and γ do.

This result corresponds in the present tunneling situation to those proven in Ref. 20 and 21 for transitions through electronic avoided crossings in the Born-Oppenheimer approximation and in Ref. 22 for wave functions obtained in semiclassical above barrier reflection situations.

Unfortunately, the explicit formulas for the emerging wave function, given by this remarkable analysis are not necessarily of practical use. By combining the analytical result with the numerical techniques, we obtain a rigorously founded spawning technique that specifies when and how to add a new family of semi-classical wavepackets during the tunneling in one dimension.

III. SPAWNING THE GAUSSIAN

We start with a linear combination (8) of wave-packets corresponding to a family of parameters q, p, Q, P and take the number of basis functions K large enough to accurately propagate the wave-function even somewhat on the far side of the barrier. The energy balance (numerical conservation of the energy) can indicate whether K should be increased. The position parameter q indicates whether the propagation time is sufficient. The theory says that the part of the wave-function beyond the barrier is approximately a Gaussian after sufficient time. Using only the parameters of the propagated wave-function q, p, Q, P, S, (c k ) k≤K we want to create the new set of parameters a, b, A, B, d 0 so that the tunneled wave-packet is approximately

d 0 ϕ 0 [a, b, A, B].
The main observation is that the coefficients c k for large indices k are responsible for the correct representation of the tunneled part of the wave-function. After tunneling has occurred, we decompose the wave-function into two parts

u(t) = v(t) + w(t) = k<K0 c k (t)ϕ k + K k=K0 c k (t)ϕ k
where w is expected to be approximately a Gaussian:

w(t, x) ≈ d 0 ϕ 0 [a, b, A, B](x) . Of course, w 2 = K k=K0 |c k | 2 = 1, but ϕ 0 [a, b, A, B]
is normed, so we may set d 0 = w . Note that the choice of K 0 is not critical, as long as the two parts of the wavefunction are well separated (i.e., the coefficients c k in the middle range of k are nearly zero). The estimated position and momentum are

a = w, xw w 2 , b = w, -iε∇w w 2 .
Using (7) and the orthogonality of the basis functions, we get

a = q + √ 2ε w 2 Re Q K k=K0+1 c k c k-1 √ k , (14) and b
= p + √ 2ε w 2 Re P K k=K0+1 c k c k-1 √ k . (15)
We use Ref. 2, equation (3.28)

ε 2 |A| 2 = w, (x -a) 2 w w 2 ,
together with the last two equalities in order to get

|A| 2 = - 2 ε (q -a) 2 + 1 w 2 |Q| 2 θ 1 + 2Re (Q 2 θ 2 ) ,
where

θ 1 = K k=K0 |c k | 2 (2k + 1) θ 2 = K-2 k=K0 c k+2 c k (k + 1)(k + 2) .
Analogously we obtain a formula for |B| 2 :

|B| 2 = - 2 ε (p -b) 2 + 1 w 2 |P | 2 θ 1 + 2Re (P 2 θ 2 ) .
At the expense of multiplying d 0 by a phase, we can choose a real A > 0. Then the second compatibility relation (4) ensures that Im B = 1/A and hence

A = |A| 2 (16)
B = A 2 |B| 2 -1 + i A . ( 17
)
The equations ( 14), ( 15), ( 16), (17) give the parameters of the spawned Gaussian in the tunneling problem. Note that the same algorithm will provide the parameters of the Gaussian that emerges in the reflection problem.

Let us point out here the weak point of the method: for an accurate spawning we may need a large number K of basis functions in the case of large initial momentum. The computing time prior to spawning in this case may be substantial. However, spawning reduces dramatically the number of basis functions that are needed after tunneling has occurred, allowing for long time propagations.

IV. NUMERICAL RESULTS

The first question to address is whether our propagation scheme can reproduce the results obtained by the TDDVR 10 . First, we have to describe the problem in Ref. 10 in terms of our parameters and of atomic units that are suitable for numerical computations. The Eckart potential is

V (x) = V 0 cosh 2 (ax) , with V 0 = 100 kJ/mol = 100 • 3.8008 • 10 -4 Eh and a = 2 Å-1 = 2 • 0.52918 alu -1 .
The parameter ε 2 is the mass ratio between the electrons and the considered system mass of 1 amu, hence ε = 0.02342. Note that this value of ε is large enough in order to address the one dimensional problem by the Fourier method on a large computational domain. This will provide us a very good approximation for the exact solution of (1). The unit of time in (1) is then 1/ε = 1.0327 fs.

The propagation starts with a Gaussian u(x, 0) = ϕ ε 0 [q(0), p(0), Q(0), P (0)](x). The initial position of the wave-packet is q(0) = -4 Å = -4/0.52918 and the initial momentum is p

(0) = ε • 20 Å-1 = ε • 20 • 0.52918. Tak- ing the initial Q(0) = 1/ √ 2 • 0.4 and P (0) = i √ 2
• 0.4 lets us exactly reproduce the case treated in Ref. 10 (see Figure 3). Note the perfect conservation of the total energy of 3.1192 • 10 -2 Eh during the propagation in Figure 4. With this choice of parameters, one can easily compare the solution given by our algorithm with that provided by the split-step Fourier method of Ref. 25, see Figure 5. In order to compute an approximation of the L 2 -and L ∞ -norm (or max-norm) of the error, both solutions are evaluated on an uniform grid with 4096 points in the computational box [-9π, 9π] that is reqired by the Fourier method. Note that the wavepacket propagation happens on the whole real line, so there is no need to impose artificial boundary conditions. Moreover, smaller ǫ would require many more frequencies (and points) in the Fourier method in order to make an accurate comparison. We chose the time-step to be 0.005, and we used 512 basis functions in the wave-packet Ansatz. The right part of Figure 5 shows the same oscillations in the tunneling probability 5 we notice that the approximation deteriorates during the tunneling, and it remains constant thereafter. To avoid approximating two independently propagating waves by only one wave-packet, we spawn a new family as described in the previous section. The left side of Figure 6 shows the L 2 -and L ∞ -norm (or max-norm) of the spawing error |w(x) -d 0 ϕ 0 [a, b, A, B](x)| at different possible spawing times during and after the tunneling. The constant K 0 was chosen to be 50, so one can cope with many fewer basis-functions by introducing the second wave-packet, see Figures 1 and2. The right side of Figure 6 and Figure 2 present the corresponding simmulation results in the case when we take as an initial wave-packet u(x, 0) = ϕ ε 3 [q(t), p(t), Q(t), P (t)](x) with the same q(0), Q(0), P (0) and slightly less momentum p(0) = 0.95 • ε • 20 Å-1 , hence with the total energy 2.9264 • 10 -2 Eh.

In conclusion, the spawning technique reduces dramatically the number of basis functions used, while keeping the error resonably small.

FIG. 1 .

 1 FIG. 1. The numerical solution (for u(0) = ϕ0) at several times: upper lines display |u| 2 and lower lines show |c k |; the dotted line represents the spawned Gaussian, the colors indicate the phase of the wavefunction 15 .

FIG. 2 .

 2 FIG. 2. The numerical solution (for u(0) = ϕ3) at several times: upper lines display |u| 2 and lower lines show |c k |; the dotted line represents the spawned Gaussian, the colors indicate the phase of the wavefunction 15 .

∞0

  |u(t, x)| 2 dx as

FIG. 3 . 24 FIG. 4 .

 3244 FIG.3. Energy-resolved transmission probabilities; the numerical computed values are compared with the analytic result[START_REF] Landau | Quantum mechanics: nonrelativistic theory[END_REF][START_REF] Billing | Advanced molecular dynamics and chemical kinetics[END_REF]