Thus, the wave packet propagation (WPP) method is success- fully applied to time dependent and time-independent problems in gas-surface interactions, ...

0 downloads 12 Views 114KB Size

No documents

arXiv:physics/0110077v1 [physics.comp-ph] 26 Oct 2001

Andrei G. BORISOV

a

and Sergei V. SHABANOV

b

a

Laboratoire des Collisions Atomiques et Mol´eculaires, Unit´e mixte de recherche CNRS-Universit´e Paris-Sud UMR 8625, Bˆ atiment 351, Universit´e Paris-Sud, 91405 Orsay CEDEX, France b Institute for Fundamental Theory, Departments of Physics and Mathematics, University of Florida, Gainesville, FL 23611, USA Abstract It is demonstrated that the wavelets can be used to considerably speed up simulations of the wave packet propagation in multiscale systems. Extremely high efficiency is obtained in the representation of both bound and continuum states. The new method is compared with the fast Fourier algorithm. Depending on ratios of typical scales of a quantum system in question, the wavelet method appears to be faster by a few orders of magnitude.

1. Owing to the fast development of the computational tools the direct solution of the time-dependent Schroedinger equation has become one of the basic approaches to study the evolution of quantum systems. Thus, the wave packet propagation (WPP) method is successfully applied to time dependent and time-independent problems in gas-surface interactions, molecular and atomic physics, and quantum chemistry [1, 2, 3, 4]. One of the main issues of the numerical approaches to the time-dependent Schroedinger equation is the representation of the wave function |Ψ(t)i of the system. Development of the pseudospectral global grid representation approaches yielded a very efficient way to tackle this problem. The discrete variable representation (DVR) [5] and the Fourier grid Hamiltonian method (FGH) [6, 7] have been widely used in time-dependent molecular dynamics [1, 2, 8], S-matrix [9, 10] or eigenvalue [11] calculations. The FGH method based on the Fast Fourier Transform (FFT) algorithm is especially popular. This is because for the mesh of N points the evaluation of the kinetic energy operator requires only NlogN operations and can be easily implemented numerically. In the standard FGH method the wave function is represented on a grid of equidistant point in coordinate and momentum space. It is well suited when the the local de Broglie wavelength does not change much over the physical region spanned by the wave function of the system [11]. There are, however, lot of examples where the system under consideration spans the physical regions with very different properties. Consider, for example, the scattering of a slow particle on a narrow and deep potential well. Despite the short wave lengths occur only in the well, in the FGH method they would determine the lattice mesh over entire physical space leading to high computational cost. In fact, pseudospectral global grid representation approaches are difficult to use in multiscale problems. This is why much of work has been devoted recently to the development of the mapping procedures in order to enhance sampling efficiency in the regions of the rapid variation of the wave function [8, 11, 12, 13, 14, 15]. Though mapping procedure, based on the variable change x = f (ξ) is very efficient in 1D, it is far from being universal. One obvious drawback 1

is that it is difficult to implement in higher dimensions. In the case of several variables, there is no simple procedure to define the mapping functions f so that the lattice would be fine only in some designated regions of space. Next, the topology of the new coordinate surfaces can be different from that of Cartesian planes. Therefore the Jacobian may vanish at some points (e.g., spherical coordinates). This leads to singularity in the kinetic energy and imposes small time step in simulations [16]. In this letter we propose a novel approach to the wave packet propagation in multiscale systems. It is based on the use of wavelets as basis functions in the projected Hilbert space. In contrast to the plane waves, the wavelets are localized in both real and momentum spaces. This characteristic property of wavelets allows us to accurately describe short wave length components of the wave function in designated spatial regions, while keeping only few basis elements in the bulk. This leads to a drastic reduction of the projected Hilbert space dimension without any loss of accuracy. In some cases it may even become possible to diagonalize the Hamiltonian matrix so that the time evolution can be elementary followed in the basis of eigenstates. We illustrate our approach by simulations of a simple multiscale one-dimensional scattering problem. A generalization to higher dimensions is straightforward and based on the standard mathematical construction of multidimensional wavelets [17]. Wavelet bases have been successfully used for the systematic treatment of the multiple length scales in the electronic structure of matter [18, 19, 20, 21] (for a comprehensive review see [22] and references therein). Here we apply, for the first time, wavelets to study the time evolution in multiscale quantum systems and demonstrate the advantages of the wavelet method over the conventional FGH method. 2. To illustrate the efficiency of the wavelet method in the wave packet propagation, we consider the 1D scattering of an electron on a narrow and deep potential well. Despite we use an electron as a projectile, procedure described below readily applies to molecular or atomic wave packets. The potential has the form Vˆ = −64 exp[−(100x)2 ] (atomic units are used throughout the paper) and supports single bound state with the energy −0.63200. We are interested in the transmission and reflection coefficients for the energies of the scattered electron within 0.2 − 1.1 range. This corresponds to a typical wave length λ = 2π which is much larger than the potential well width. To calculate the scattering properties of the system we use the Gaussan wave packet impinging on a potential well. The calculation requires a typical size of the box x ∈ [−60, 60], where 18 a.u. from each side are allocated for an optical potential Vˆopt that absorbs reflected and transmitted waves [7]. We first introduce a uniform lattice and simulate the wave packet propagation by means of the FGH method. Split-operator technique [6, 16] is used to calculate the action of the ˆ = Vˆ + Vˆopt ): evolution operator (W ˆ + Vˆopt ) |Ψ(t)i |Ψ(t + ∆t)i = exp(−i∆t H h

i

ˆ /2) exp(−i∆tTˆ) exp(−i∆tW ˆ /2) |Ψ(t)i + O(∆t3 ) = exp(−i∆tW

(1)

Numerical convergence is obtained with N = 214 points at the grid and a time step ∆t = 0.0002. Calculated transmission and reflection coefficients are shown in Fig. 1. The time 2

evolution of the wave packet is presented in Fig.2. The colors represent the magnitude of the wave packet in the logarithmic scale: −10 ≥ ln |Ψ(t, x)| ≥ 0.3. As the color changes from red to violet, ln |Ψ(t, x)| decreases from 0.3 to −10. Without the potential present, we would see a colored ray (a trajectory of a free particle is a straight line) spreading for larger values of t. To compare the above results with our new wavelet method, we take the Daubechies wavelets D10 which have ten first vanishing moments and the filter of length 20 [17]. The basis is generated by the scaling function φ(x) whose support lies in the interval x ∈ [0, 19]. The lowest resolution level is set so that the initial wave packet is reproduced with high √ accuracy in the orthonormal basis of the functions φ1,j (x) = 2φ(2x √ − j) where j runs over integers inside the interval (−120, 120]. The wavelets ψk,j (x) = 2k ψ(2k x − j) are used to describe short wave length components of the wave function in the vicinity of the potential well. For a moment, they should be regarded as a special orthonormal basis with compact support and the following properties. If the basis of scaling functions φ1,j (x) spans well functions whose Fourier transforms are concentrated in a wave length band around λ1 , the corresponding wavelets from the kth resolution level form a good basis for a band centered at 2−k λ1 . In our case we have used k = 1, 2, 3, 4. The number of wavelets needed on each resolution level is determined by the potential width (by its shape, in general). The technical details are explained in Section 3. It appears necessary to take 20 wavelets on each resolution level. Thus, all together we use only 240 + 4 · 20 = 320 coefficients (or, equally basis functions) to describe the wave packet propagation instead of 214 in the fast Fourier method. Due to such a tremendous reduction of the size of the problem, the Hamiltonian matrix can be directly diagonalized and the time dependent wave function can be easily obtained as: ˆ

|Ψ(t + ∆t)i = e−iVopt ∆t

P

E

e−iE∆t |Ei hE|Ψ(t)i ,

(2)

ˆ

(3)

where |Ei and E stands for the eigenvector and eigenvalue of the Hamiltonian, respectively. Convergent results are obtained with time step dt = 0.025, i.e., much larger than in the fast Fourier method above. The results are given on Figs. 1 and 3. There is a perfect coincidence of the fast Fourier and our wavelet results, while the wavelet method is nearly 100 times faster. For the sake of comparison, we also used Eq. (1) to simulate the time evolution. The same accuracy is achieved for a time step comparable to the time step in FGH method. However there is an alternative, much better splitting scheme thanks to the localization properties of ˆ =H ˆ1 + H ˆ 2 where H ˆ 1 contains matrix elements of the wavelets. Consider the decomposition H ˆ 2 corresponds to the free basis elements localized in the vicinity of the potential well, while H space and, therefore, contains matrix elements only of the scaling functions. The propagation can then be done accordingly: ˆ

ˆ

|Ψ(t + ∆t)i = e−iVopt ∆t e−iH2 ∆t/2

nP

−iE1 ∆t E1 e

o

|E1 i hE1 | e−iH2 ∆t/2 |Ψ(t)i ,

ˆ 1 . In the spit method where |E1 i and E1 stands for the eigenvector and eigenvalue of the H (1) the error depends on the spectral range of operators involved. For larger spectral ranges, the errors get larger. The conventional way to cope with this problem is to reduce the time step. The approach (3), which we call the “wavelet tower diagonalization”, offers a better 3

ˆ 1, H ˆ 2 ] is small and alternative. The convergence here is drastically improved because (a) [H ˆ (b) the operator H1 with a large spectral range can be diagonalized. Both the properties are hardly achievable without the wavelet basis. In our example the convergence is reached for a time step dt = 0.01 (vs dt = 0.0002 in (1)). Corresponding results are presented in Fig. 1. This approach applies even better to evolving heavy-particle wave packets because the Chebyshev [16, 23] or Lanczos [16, 24] schemes can be used. These schemes are, first, more efficient than the split method and, second, allow one to take full advantage of the sparse structure of the Hamiltonian in the wavelet basis (see below). Note also that by making the potential depth greater so that the minimal wave length becomes twice shorter, we would have to increase the lattice size by factor two in the Fourier method, thus increasing the number of coefficients needed to describe the wave packet by 214 (!), while in our wavelet approach one more resolution level is to be added, implying only 20 extra wavelet coefficients in the wave packet decomposition. Thus, the wavelet approach offers a systematic and easy way to improve the accuracy of simulations. 3. The description of technical details of our approach is limited to essential practical steps necessary to reproduce our results and to apply the algorithm to new systems. A general theoretical analysis of wavelets bases in multidimensional spaces can be found in [17]. So we only discuss the Hilbert space L2 of square integrable functions of a real variable x, R dx|Ψ(x)|2 < ∞. P (i). The scaling function φ(x) is defined by the equation φ(x) = 2 l hl φ(2x − l) (called the scaling relation), where the coefficients hl are called a filter. For a finite filter, the scaling function has compact support. Define φn,j (x) = 2n/2 φ(2n x − j) for all integers j and n. The filter hl satisfies an equation obtained by combining the scaling relation with the required orthonormality condition of the scaling functions φn,j for fixed n. Given a filter, numerical values of φ(x) can be generated by an iteration procedure [17, 25]. (ii). The subspace of L2 spanned by φn,j is denoted Vn . An important property Vn is that Vn−1 ⊂ Vn and the projection of any function Ψ from L2 onto Vn converges to Ψ in the L2 sense as n → ∞. Consider an orthogonal decomposition: Vn+1 = Vn ⊕ Wn . There exists a function ψ ∈ W0 called the mother wavelet such that ψn,j (x) = 2n/2 ψ(2n x − j) form an orthonormal basis in Wn . The numerical values of ψ can be generated from the scaling relation P ψ(x) = 2 l (−1)l h1−l φ(2x − l). A finite dimensional subspace of Vn is used to approximate L2 in simulations. By construction, the functions φl,j , l < n, and ψk,j for k = l, l + 1, ..., n − 1 form an orthonormal basis in Vn : hψk,j |ψk′ ,i i = δkk′ δij , hφl,j |φl,i i = δij , and hψk,j |φl,i i = 0. (iii). From the definition of φl,j and ψk,j it is clear that the index j indicates the position of support of the scaling function or wavelet. The index k of ψk,j can be understood as follows. Let the Fourier transform of the mother wavelet ψ be peaked at a momentum p0 . Then the Fourier transform of ψk,j is peaked at the momentum pk = 2k p0 . Since the wavelets with different k are orthogonal, the coefficients dk,j = hψk,j |Ψi determine relative amplitudes of successively shorter wave length components of Ψ in the vicinity of x = 2−k j as k increases. For this reason, the index j is called a position index, and k is called a resolution level. (iv). From the physical properties of a system one can estimate a wave length band, λ ∈ [λmin , λmax ], required for simulations. The lowest resolution level l is identified by the condition that φl,j span functions with Fourier components in the vicinity of λmax . The 4

necessary maximal resolution level is determined by n ≈ log2 (λmax /λmin ). In our example l = 1 and n = 5. If N scaling functions are needed (to cover a given physical volume), then, in general, on every wavelet resolution level there will be Jk = 2k−l N basis functions. The wave packet is decomposed in the corresponding basis of Vn Ψ(t) =

N X

s1,j (t)φ1,j +

Jk n−1 XX

dk,j (t)ψk,j .

(4)

k=1 j=1

j=1

The decomposition coefficients as functions of time are to be found from the Schroedinger equation. Yet, the total number of coefficients in Ψ(t) is about the same as that in the uniform finite lattice approach. The advantage of using wavelet bases becomes significant whenever the volume of regions where short wave lengths can appear during the time evolution is much smaller than the total physical volume of the system. This allows one to substantially reduce Jk in Eq.(4) by taking higher resolution level wavelets only where they are needed. In our case we need only 20 wavelets at each resolution level. We shall refer to all the wavelets needed in the vicinity of one local minima of the potential as a wavelet tower. Higher tower floors correspond to higher resolution levels k. Each wavelet is thought of as a building block of width equal 2−k . (v). An important part of the algorithm is the projection of the Hamiltonian onto the wavelet towers. Matrix elements of the potential as well as the initial coefficients sl,j = hφl,j |Ψi, dk,j = hψk,j |Ψi are computed via Riemann sums. The matrix elements of the kinetic energy are computed using the the fast Fourier transform to obtain the second derivative of the function. The Hamiltonian matrix is sparse because of compact support of the basis functions. For instance, in our example hφn,j |H|φn,ii = 0, if |i − j| ≥ 20. Finally, the time evolution can be computed by standard techniques such as split, Lanczos, or Chebychev methods [16, 23, 24]. If after the projection on wavelet towers, the Hamiltonian matrix is small enough for a direct diagonalization, Eq. (2) can be used for the propagation. In conclusion, we have demonstrated that wavelet bases can be extremely efficient for solving the time-dependent Schroedinger equation for multiscale systems. The bound and continuum states are accurately described with a much less number of basis functions as would be required in the standard Fourier-grid methods. In our example we were able to reduce the basis size by a factor of 50, which allowed us to scale down the computation time by a factor of 100. The method can easily be implemented in higher dimensions where a wavelet basis is built by taking the direct product of one-dimensional wavelets [17]. The wavelet approach becomes more efficient as ratios of typical scales get larger. Finally we would like to stress that the wavelet method is ideally suitable for simulating wave packet propagation in multiscale problems with complicated form of the potential. This is because the wavelet towers can be custom designed for any topology of the potential minima (e.g., potential valleys) regardless of the dimensionality of the system. Acknowledgments. S.V.S. thanks the LCAM of the University of Paris-Sud for warm hospitality. We are also grateful to Victor Sidis and Jean Pierre Gauyacq for support of this project, and to John Klauder for reading the manuscript and useful comments. 5

References [1] G.-J. Kroes, Progress in Surf. Sci. 60, 1 (1999) [2] M.H. Beck, A. J¨ackle, G.A. Worth, H.-D. Meyer, Phys. Reports 324, 1 (2000) [3] J.Z.H. Zhang, Theory and Application of Quantum Molecular Dynamics, (World Scientific, New Jersey, 1999) [4] R.N. Zare, Science 279, 1875 (1998) [5] J.V. Lill, G.A. Parker, J.C. Light, Chem. Phys. Lett. 89, 483 (1982) [6] M.D. Feit, J.A. Fleck, Jr., and A. Steiger, J. Comput. Phys. 47, 418 (1982) [7] D. Kosloff and R. Kosloff, J. Comput. Phys. 52, 35 (1983) [8] R. Kosloff in Time Dependent Quantum Molecular Dynamics, NATO ASI series, Series B: Physics, edited by J. Broeckhove and L. Lathouwers (Plenum, NY 1992), Vol 299, p. 97 [9] N. Rom, J.W. Pang, and D. Neuhauser, J. Chem. Phys. 105, 10436 (1996) [10] V.A. Mandelshtam and H.S. Taylor, J. Chem. Phys. 102, 7390 (1995) [11] E. Fattal,R. Baer and R. Kosloff, Phys. Rev. E 53, 1217 (1996) [12] D. Lemoine, Chem. Phys. Lett. 320, 492 (2000) [13] F. Gygi, Phys. Rev. B 48, 11692 (1993) [14] I. Tuvi and Y.B. Band, J. Chem. Phys. 107, 9079 (1997) [15] V. Kokoouline, O. Dulieu, R. Kosloff, F. Masnou-Seeuws, J. Chem. Phys. 110, 9865 (1999) [16] C. Leforestier et al, J. Comput. Phys. 94, 59 (1991) [17] I. Daubechies, Ten Lectures on Wavelets, (SIAM, Philadelphia, 1992) [18] K. Cho, T.A. Arias, J.D. Joannopoulos, P.K. Lam, Phys. Rev. Lett. 71, 1808 (1993) [19] J.P. Modisette, P. Nordlander, J.L. Kinsey, and B.R. Johnson, Chem. Phys. Lett. 250, 485 (1996) [20] S. Wei and M.Y. Chou, Phys. Rev. Lett. 76, 2650 (1996) [21] C.J. Tymczak and X.-Q. Wang, Phys. Rev. Lett. 78, 3654 (1997) [22] T. A. Arias, Rev. Mod. Phys. 71, 267 (1999)

6

[23] H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984) [24] T.J. Park and J.C. Light, J. Chem. Phys. 85, 5870 (1986) [25] W.H. Press et al, Numerical Recipies in Fortran, Second Edition, (Cambridge University Press, Cambridge, 1994), p. 584

Figure captions Fig. 1. Calculated transmission (black) and reflection (red) coefficients. Circles: The standard approach based on the split propagation and Fourier grid with uniform mesh. Solid curves: The wavelet approach with Hamiltonian matrix diagonaliztion. Triangles: The split propagation with the “wavelet tower diagonalization” as in Eq. (3).

7

transmission and reflection coefficients

0.8

0.7

0.6

0.5

0.4

0.3

0.2 0.2 0.3 0.4 0.5

0.6 0.7 0.8

energy (a.u.)

0.9 1.0

1.1

This figure "f2.gif" is available in "gif" format from: http://arxiv.org/ps/physics/0110077v1

This figure "f3.gif" is available in "gif" format from: http://arxiv.org/ps/physics/0110077v1