###### 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.

The wave packet propagation using wavelets

Andrei G. BORISOV and Sergei V. SHABANOV

Laboratoire des Collisions Atomiques et Moléculaires,

Unité mixte de recherche CNRS-Université Paris-Sud UMR 8625,

Bâtiment 351, Université Paris-Sud, 91405 Orsay CEDEX, France

Institute for Fundamental Theory, Departments of Physics

and Mathematics, University of Florida, Gainesville, FL 23611, USA

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 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 points the evaluation of the kinetic energy operator requires only 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 is very efficient in 1D, it is far from being universal. One obvious drawback 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 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 (atomic units are used throughout the paper) and supports single bound state with the energy . We are interested in the transmission and reflection coefficients for the energies of the scattered electron within range. This corresponds to a typical wave length 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 , where 18 a.u. from each side are allocated for an optical potential 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 evolution operator :

(1) | |||||

Numerical convergence is obtained with points at the grid and a time step . Calculated transmission and reflection coefficients are shown in Fig. 1. The time evolution of the wave packet is presented in Fig.2. The colors represent the magnitude of the wave packet in the logarithmic scale: . As the color changes from red to violet, decreases from to . 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 .

To compare the above results with our new wavelet method, we take the Daubechies wavelets which have ten first vanishing moments and the filter of length 20 [17]. The basis is generated by the scaling function whose support lies in the interval . The lowest resolution level is set so that the initial wave packet is reproduced with high accuracy in the orthonormal basis of the functions where runs over integers inside the interval . The wavelets 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 spans well functions whose Fourier transforms are concentrated in a wave length band around , the corresponding wavelets from the th resolution level form a good basis for a band centered at . In our case we have used . 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 coefficients (or, equally basis functions) to describe the wave packet propagation instead of 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:

(2) |

where and stands for the eigenvector and eigenvalue of the Hamiltonian, respectively. Convergent results are obtained with time step , 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 wavelets. Consider the decomposition where contains matrix elements of the basis elements localized in the vicinity of the potential well, while corresponds to the free space and, therefore, contains matrix elements only of the scaling functions. The propagation can then be done accordingly:

(3) |

where and stands for the eigenvector and eigenvalue of the . In the spit method (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 alternative. The convergence here is drastically improved because (a) is small and (b) the operator 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 (vs 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 (!), 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 of square integrable functions of a real variable , .

(i). The scaling function is defined by the equation (called the scaling relation), where the coefficients are called a filter. For a finite filter, the scaling function has compact support. Define for all integers and . The filter satisfies an equation obtained by combining the scaling relation with the required orthonormality condition of the scaling functions for fixed . Given a filter, numerical values of can be generated by an iteration procedure [17, 25].

(ii). The subspace of spanned by is denoted . An important property is that and the projection of any function from onto converges to in the sense as . Consider an orthogonal decomposition: . There exists a function called the mother wavelet such that form an orthonormal basis in . The numerical values of can be generated from the scaling relation . A finite dimensional subspace of is used to approximate in simulations. By construction, the functions , , and for form an orthonormal basis in : , , and .

(iii). From the definition of and it is clear that the index indicates the position of support of the scaling function or wavelet. The index of can be understood as follows. Let the Fourier transform of the mother wavelet be peaked at a momentum . Then the Fourier transform of is peaked at the momentum . Since the wavelets with different are orthogonal, the coefficients determine relative amplitudes of successively shorter wave length components of in the vicinity of as increases. For this reason, the index is called a position index, and is called a resolution level.

(iv). From the physical properties of a system one can estimate a wave length band, , required for simulations. The lowest resolution level is identified by the condition that span functions with Fourier components in the vicinity of . The necessary maximal resolution level is determined by . In our example and . If scaling functions are needed (to cover a given physical volume), then, in general, on every wavelet resolution level there will be basis functions. The wave packet is decomposed in the corresponding basis of

(4) |

The decomposition coefficients as functions of time are to be found from the Schroedinger equation. Yet, the total number of coefficients in 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 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 . Each wavelet is thought of as a building block of width equal .

(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 , 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 , if . 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.

## References

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).