Recently we introduced a phase space approach for solving the time-independent Schrödinger equation using a periodic von Neumann basis with bi-orthogonal exchange (pvb) [A. Shimshovitz and D. J. Tannor, Phys. Rev. Lett.109, 070402 (2012)

]. Here we extend the approach to allow a wavelet scaling of the phase space Gaussians. The new basis set, which we call the wavelet pvb basis, is simple to implement and provides an appealing alternative to other wavelet approaches. For the 1D Coulomb problems tested in this paper, the method reduces the size of the basis relative to the Fourier grid method by a factor of 13–60. The savings in basis set size is predicted to grow steeply as the dimensionality increases.

The term wavelet refers to a broad class of basis sets in which the basis functions are all replicas of one “mother” function with different characteristic scales. In contrast to Fourier analysis, the wavelet functions are localized both in time and frequency and therefore are a natural tool for the analysis of signals where combined time-frequency resolution is required.1,2 In addition, wavelets are useful for signal compression since they provide a compact representation of signals with different frequency components at different times.3–6 For the same reason, wavelets can provide an efficient representation of wavefunctions in quantum mechanics if the wavefunction varies rapidly in some places and slowly in others.7–11 

To some extent, wavelets may be viewed as the generalization of a phase space localized basis set introduced by von Neumann (vN) in 193112 and independently by Gabor in 1946 in the context of signal processing.13 In the vN basis, the phase space localization is achieved by shifting a fiducial Gaussian in coordinate and momentum space, resulting in basis functions that are centered on equidistant unit cells in phase space. Because in practice the phase space lattice must always be truncated, there have been longstanding convergence problems with the vN basis. We have recently shown that two modifications are needed to achieve high accuracy and efficiency with the von Neumann basis: (1) to use a periodic version of the vN basis (pvN) and (2) to exchange the role of the pvN basis with its bi-orthogonal basis, so that localized von Neumann functions are used to calculate the coefficients, rather than serving as the basis14 (periodic von Neumann with biorthogonal exchange or pvb). In the present paper we combine the pvb method with the insight that the division into equidistant unit cells and equal width Gaussians is not necessary. In particular, we focus on a 1-parameter scaling transformation that links the width and the spacing of the Gaussians and show that by adding this freedom we can significantly improve the efficiency of the calculation. We call the new basis wavelet pvb or wpvb.

The wpvb basis set has several significant advantages over other wavelet methods that have been used in the quantum context: (1) Simplicity: The construction of the Hamiltonian elements is done simply, accurately, and efficiently using matrix multiplication, while in conventional wavelet methods one needs sophisticated algorithms to get approximations to the Hamiltonian matrix elements. (2) Phase space interpretation: Because the wpvb method is based on Gaussians, the location of each basis function in phase space is well-defined. (3) Locality: The Gaussian function used in the wpvb method is optimally concentrated in both momentum and coordinate and therefore expected to be the most suitable for a basis set with high phase space resolution. (4) Flexibility: The wpvb basis does not require a rigid scaling relationship, and in principle can be tailored to fit arbitrarily complicated phase space geometries. Conventional wavelet methods do not have this flexibility. (5) Precision: The wpvb is exact for functions that are band limited and of finite support and converges exponentially for wavefunctions that are semilocal in phase space. We present results below for which the convergence in the wpvb method is competitive and even slightly faster than Daubechies's wavelets2 (Fig. 3(b)).

Our starting point is the Fourier grid (FG) representation15–17 in which the wave function ψ(x) is sampled on some finite, evenly spaced grid {xi} such that

$\psi (x)\break\rightarrow \tilde{\psi }(x_i)$
ψ(x)ψ̃(xi)⁠, xi = (i − 1)Δ, i = 1, …, N, Δ is the spacing between grid points and L = NΔ is interval spanned by the grid. This representation is exact in the interval [0, L] for any wavepacket ψ(x) with support between [0, L] and band limited between ( − Pmax, Pmax) in momentum, where
$P_{\rm {max}}=\frac{\pi \hbar }{\Delta }$
P max =πΔ
. The reconstruction of ψ(x) from its samples is given by

\begin{eqnarray}\psi (x)=\sum _{n=1}^{N} \psi (x_n)\theta _{n}(x),\end{eqnarray}
ψ(x)=n=1Nψ(xn)θn(x),
(1)

where θn(x) is the Dirichlet or periodic sinc function.18,19 For wavefunctions that are semilocalized in phase space, the FG representation converges exponentially with the number of grid points.15 

The pvb basis is formally fully equivalent to but more flexible than the FG basis. It is based on the von Neumann functions,12 

\begin{eqnarray}g_{nl}(x)=\left(\frac{2}{a^2}\right)^\frac{1}{4}\exp \left(-\frac{\pi }{a^2}(x-na)^{2}+l\frac{2\pi i}{a}(x-na)\right),\nonumber\\\end{eqnarray}
gnl(x)=2a214expπa2(xna)2+l2πia(xna),
(2)

where the indices n and l take a finite number of values: n = 0, …, Nx − 1 and

$l=\frac{-N_p}{2}+1$
l=Np2+1⁠, …,
$\frac{N_p}{2}$
Np2
. This produces N Gaussian basis functions centered at
$(na,\frac{lh}{a})$
(na,lha)
with spacing a and
$\frac{h}{a}$
ha
in x and p, respectively (see Fig. 1(a)). In order to obtain equivalence with the FG representation, we use the set {gm(x)}, m = 1, …, N, to produce a linear combination of the Dirichlet functions as follows:

\begin{eqnarray}\tilde{g}_m(x)=\sum _{i = 1}^{N} \theta _i(x)g_m(x_i).\end{eqnarray}
g̃m(x)=i=1Nθi(x)gm(xi).
(3)

Since the {θi(x)} are periodic in L, so are the

$\lbrace \tilde{g}_m(x)\rbrace$
{g̃m(x)}⁠, hence we call this new basis the periodic vN basis. The transformation equation (3) is invertible, so the pvN basis set spans the same Hilbert space as the FG representation. However, the basis
$\lbrace \tilde{g}_m(x)\rbrace$
{g̃m(x)}
is non-orthogonal so the expansion of ψ(x) (Eq. (1)) takes the following form:

\begin{equation}\vert \psi \rangle =\sum _{n = 1}^{N} \sum _{m = 1}^{N} \vert \tilde{g}_m \rangle (S^{-1})_{mn} \langle \tilde{g}_n\vert \psi \rangle,\end{equation}
|ψ=n=1Nm=1N|g̃m(S1)mng̃n|ψ,
(4)

where

\begin{eqnarray}\langle \tilde{g}_n\vert \psi \rangle &=& \sum _{n=1}^Ng^{*}_i(x_n)\psi (x_n),\end{eqnarray}
g̃n|ψ=n=1Ngi*(xn)ψ(xn),
(5)

S is the N × N overlap matrix with elements

\begin{eqnarray}S_{ij}= \langle \tilde{g}_i \vert \tilde{g}_j\rangle =\sum _{n=1}^Ng^{*}_i(x_n)g_j(x_n),\end{eqnarray}
Sij=g̃i|g̃j=n=1Ngi*(xn)gj(xn),
(6)

and we have used the orthonormality of the {θ} functions. Note that only the {gi} (not

$\lbrace \tilde{g}_i\rbrace$
{g̃i}⁠) enter into Eqs. (5) and (6). Although Eqs. (5) and (6) can be viewed as discrete approximations to continuous integrals, paradoxically this “approximation” gives orders of magnitude improvement in accuracy over calculating the integrals exactly.14 

Since the pvN functions are localized in phase space (Fig. 1), we expect that the factor

$\langle \tilde{g}_n\vert \psi \rangle$
g̃n|ψ in Eq. (4) will be nearly zero for some of the
$\tilde{g}_n$
g̃n
functions. We can take advantage of this by defining a basis {bi} bi-orthogonal to the periodic von Neumann basis (henceforth “pvb”) as

\begin{eqnarray}\vert \tilde{b}_i\rangle =\sum _{j = 1}^{N} \vert \tilde{g}_j\rangle (S^{-1})_{ji}.\end{eqnarray}
|b̃i=j=1N|g̃j(S1)ji.
(7)

Comparing Eq. (7) with Eq. (4), we find that

\begin{eqnarray}\vert \psi \rangle =\sum _{n = 1}^{N} \vert \tilde{b}_n\rangle \langle \tilde{g}_n\vert \psi \rangle \equiv \sum _{n = 1}^{N} \vert \tilde{b}_n\rangle c_n.\end{eqnarray}
|ψ=n=1N|b̃ng̃n|ψn=1N|b̃ncn.
(8)

If M of the overlaps

$\langle \tilde{g}_n\vert \psi \rangle$
g̃n|ψ are zero, M of the cn coefficients will be zero and the corresponding basis functions can be discarded. Note the counterintuitive property that
$\lbrace \tilde{b_i}\rbrace$
{bĩ}
functions may be discarded even if they are delocalized if their coefficients
$\langle \tilde{g}_n\vert \psi \rangle$
g̃n|ψ
are nearly zero.

FIG. 1.

(a) N = 16 vN unit cells. The vN basis functions are Gaussians located at the center of each unit cell. (b) wpvb basis set. The Gaussian basis functions may have different widths in x and p retaining the unit cell area equal to h.

FIG. 1.

(a) N = 16 vN unit cells. The vN basis functions are Gaussians located at the center of each unit cell. (b) wpvb basis set. The Gaussian basis functions may have different widths in x and p retaining the unit cell area equal to h.

Close modal

An important application of the pvb basis set is to eigenvalue problems. In the Fourier Grid Hamiltonian (FGH)/sinc-DVR method,16,17,20 one constructs a Hamiltonian matrix

$\rm {\textbf {H}}^{\rm FGH}$
H FGH in which the potential matrix is diagonal with elements Vij = V(xiij and the kinetic energy matrix is given by

\begin{eqnarray}T^{\rm FGH}_{ij} =\frac{\hbar ^{2}}{2M} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{K^{2}}{3}(1+\frac{2}{N^{2}}), & \mbox{if} \quad i = j, \\[6pt]\frac{2K^{2}}{N^{2}}\frac{(-1)^{j-i}}{\rm {sin}^{2}(\pi \frac{j-i}{N})} & \mbox{if} \quad i \ne j. \end{array}\right.\end{eqnarray}
Tij FGH =22MK23(1+2N2),ifi=j,2K2N2(1)ji sin 2(πjiN)ifij.
(9)

We can transform HFGH to the pvb representation by

\begin{eqnarray}H^{\rm pvb}_{ij}=\sum _{m = 1}^{N}\sum _{n = 1}^{N} b_i^{*}(x_m)H^{\rm FGH}_{mn}b_j(x_n),\end{eqnarray}
Hij pvb =m=1Nn=1Nbi*(xm)Hmn FGH bj(xn),
(10)

or

$\rm {\textbf {H}}^{pvb}=\rm {\textbf {B}}^\dagger \rm {\textbf {H}}^{\rm {FGH}}\rm {\textbf {B}}$
H pvb =BH FGH B⁠. Again, only the {bi} (not the
$\lbrace \tilde{b}_i\rbrace$
{b̃i}
) enter into Eq. (10). Note that the transformation matrix B is non-unitary; hence we have to solve the generalized eigenvalue problem
$\rm {\textbf {H}}^{\rm pvb}\rm {\textbf {U}}=\rm {\textbf {S}}^{\rm pvb}\rm {\textbf {UE}}$
H pvb U=S pvb UE
, where
$\rm {\textbf {S}}^{\rm pvb}$
S pvb
is the overlap matrix in the pvb basis set which can be shown to be equal to
$\rm {\textbf {S}}^{-1}$
S1
. Although we have gone from a standard to a generalized eigenvalue problem, since the pvb basis is related to the FG basis by an invertible transformation, we obtain the same eigenenergies as for the original problem. The advantage of the new formulation is that pvb functions centered outside the classically allowed phase space can be eliminated with negligible error.

So far we have used classical phase space considerations to determine which basis functions to remove but not the shape of the basis functions themselves. However, inspection of Eq. (1) shows that the formalism can be extended to any localized basis functions. Mathematically, the only requirement for equivalence with the FGH method is that the N functions be linearly independent, although in practice, to obtain numerical stability the functions should correspond to a complete tiling of the FG phase space using N tiles, each satisfying ΔxΔp = h (see, e.g., Fig. 1(b)). As a result, we can follow the same procedure as above but using a scaling relation for the Gaussians that better matches the phase space of the problem. An important application of this scaled basis is to problems with Coulomb interactions, where the classical phase space extends to large x but has sharp features in p near the origin (Fig. 2). Representation of such a geometry is extremely difficult with the FGH method. The pvb basis is more efficient but still requires many Gaussians to span the phase space. We therefore define a new basis (henceforth wavelet pvN or “wpvN”) in which the width in x of each Gaussian gets smaller by some factor b for each level in p (Fig. 1(b)). Equation (2) then turns into

\begin{eqnarray}g_{nl}(x)=\left(\frac{2\alpha_l}{\pi}\right)^\frac{1}{4}\exp \left(-\alpha_l(x-x_{nl})^{2}+\frac{ip_l}{\hbar }(x-x_{nl})\right),\nonumber\\\end{eqnarray}
gnl(x)=2αlπ14expαl(xxnl)2+ipl(xxnl),
(11)

where al = ab|l| − 1,

$x_{nl}=(n-\frac{1}{2})a_l+x_1$
xnl=(n12)al+x1⁠,
$p_l{=}\frac{2\pi \hbar }{a_l}(\frac{b^l-1}{b-1}{-}\frac{1}{2})$
pl=2πal(bl1b112)
(l > 0), pl = −pl (l < 0), and x1 is the first point of the Fourier grid. The indices n, l take the values
$n=1,\ldots,N_x^{l}$
n=1,...,Nxl
, l = ±1, …,±Np, where
$N_x^l$
Nxl
is the number of Gaussians at each p level, given by the ratio
$N_x^{l}=\frac{L}{a_l}$
Nxl=Lal
and Np is the number of scales. The total number of Gaussians is
$\sum _{l}N_x^{l}=N.$
lNxl=N.
We call the biorthogonal basis “wpvb.” We first tested the wpvb basis set on the 1D Coulomb problem V(x) = −1/|x|. This problem is quite challenging for the FGH method since the singularity at x = 0 causes the classical phase space to extend to infinite p. We calculated the third eigenvalue (
$E_3^{{\rm {exact}}}=-\frac{1}{2(3)^2}$
E3 exact =12(3)2
) using the FGH, pvb, and wpvb methods, and examined the accuracy as a function of basis set size (Fig. 3). We used a Fourier grid of 214 points between x1 = −93.4 and xn = 93.4. The pvb basis is a product of 128 Gaussians in x and 128 Gaussians in p. For the wpvb we used the scaling given in Eq. (11) with
$a=\frac{L}{32}$
a=L32
,
$b=\frac{1}{2}$
b=12
,
$\alpha_l=\frac{1}{2a_l^2}$
αl=12al2
, and Np = 8. Pictorially, we divide the positive and the negative momentum parts of the phase space into eight inequivalent cells where the width of each cell in x is halved at each step in p. In each p segment l the whole x range is divided into
$N_x^l$
Nxl
parts, keeping the relation ΔxΔp = h. To obtain 7 digits of accuracy, the FGH method required 8000 grid points, the pvb required only 460 basis functions, and the wpvb required only 130 basis functions. The efficiency of the wpvb can be understood by examining its phase space structure (Fig. 4). All the positive p cells required for 7 digits of accuracy are marked (230 pvb and 65 wpvb cells). By scaling the wpvb cells, the narrow strip in momentum around x = 0 is spanned efficiently by the corresponding basis functions.

FIG. 2.

(a) The one-dimensional Coulomb potential as a function of x. (b) the phase space contour of

$E_3=-\frac{1}{18}$
E3=118⁠.

FIG. 2.

(a) The one-dimensional Coulomb potential as a function of x. (b) the phase space contour of

$E_3=-\frac{1}{18}$
E3=118⁠.

Close modal
FIG. 4.

Phase space cells needed to span the (positive p) classical phase space for the third eigenstate of the Coulomb potential. The negative p cells (not shown) are the mirror image. The cells correspond to the (pvb or wpvb) basis functions, both in their centers and in their aspect ratios. (a) pvb cells and (b) wpvb cells.

FIG. 4.

Phase space cells needed to span the (positive p) classical phase space for the third eigenstate of the Coulomb potential. The negative p cells (not shown) are the mirror image. The cells correspond to the (pvb or wpvb) basis functions, both in their centers and in their aspect ratios. (a) pvb cells and (b) wpvb cells.

Close modal
FIG. 3.

(a) The error in E3 for the Coulomb potential as a function of basis set size. FGH (dashed), pvb (circles), and wpvb (solid). (b) The error in E6 of the double well Coulomb potential as a function of basis set size. FGH (dashed), pvb (circles), wpvb (solid), and the results of Johnson et al. using Daubechies's wavelets (triangles).

FIG. 3.

(a) The error in E3 for the Coulomb potential as a function of basis set size. FGH (dashed), pvb (circles), and wpvb (solid). (b) The error in E6 of the double well Coulomb potential as a function of basis set size. FGH (dashed), pvb (circles), wpvb (solid), and the results of Johnson et al. using Daubechies's wavelets (triangles).

Close modal

Next, we tested the wpvb on a model for

${\rm H}_2^+$
H2+ comprising two soft Coulomb potentials:

\begin{eqnarray}\hspace*{-1pc}V(x)= -\frac{1}{\sqrt{(x-1/2)^2+\alpha ^2}}-\frac{1}{\sqrt{(x+1/2)^2+\alpha ^2}},\end{eqnarray}
V(x)=1(x1/2)2+α21(x+1/2)2+α2,
(12)

where α = 0.1 and ℏ = m = 1. Although this potential has no singularities, the classical phase space has sharp peaks in p which make the pvb (and FGH) representations inefficient (Fig. 5). We calculated the sixth eigenvalue (E6 = −0.24825962) using the FGH, pvb, and wpvb methods and studied the accuracy as a function of the basis set size. We used a grid of N = 1984 points between x1 = −34.96 and xN = 34.95. The pvb contains the product of 62 Gaussians in x with 32 Gaussians in p. For the wpvb we adopted a scaling similar to the one used for the previous problem, using Eq. (11) with

$a=\frac{L}{32}$
a=L32⁠,
$b=\frac{1}{2}$
b=12
,
$\alpha_l=\frac{1}{2a_l^2}$
αl=12al2
, and Np = 5. We also compared with a similar calculation done by Johnson et al. using Daubechies's wavelets (Fig. 3(b)). To obtain 10 digits of accuracy, the FGH method required 1560 grid points, the pvb required only 210 basis functions, and the wpvb required only 120 basis functions.

FIG. 5.

(a) The double well Coulomb potential as a function of x. (b) The phase space contour of E6 = −0.24825962.

FIG. 5.

(a) The double well Coulomb potential as a function of x. (b) The phase space contour of E6 = −0.24825962.

Close modal

The efficiency of the wpvb can be understood by examining its phase space structure (Fig. 6). All the positive p cells required for 10 digits of accuracy are marked (105 pvb and 60 wpvb cells). By scaling the wpvb cells, the sharp double-peaked structure in p near the origin is spanned efficiently.

FIG. 6.

Same as Fig. 4 but for the double well Coulomb potential. (a) pvb cells and (b) wpvb cells. Note that the wpvb functions reflect the double peak structure of the phase-space.

FIG. 6.

Same as Fig. 4 but for the double well Coulomb potential. (a) pvb cells and (b) wpvb cells. Note that the wpvb functions reflect the double peak structure of the phase-space.

Close modal

It is worth comparing the wavelet approach presented here with other wavelet approaches and with the windowed Fourier transform. (1) As in the Morlet wavelet (M-wavelet) approach, in the wpvN basis the fiducial function that is shifted and scaled is a Gaussian.21 But as opposed to the M-wavelets, the wpvN basis is defined on a discrete lattice and has periodic boundary conditions. This gives exact equivalence with the FG method without oversampling. Moreover, to achieve compression we exchange the wpvN and its biorthogonal basis, wpvb, something that is not done with M-wavelets. (2) In the Daubechies wavelet (D-wavelets) approach, the basis functions are orthogonal. Orthogonalization produces a compromise between localizing the basis functions and localizing the coefficients. Since the wpvN is a non-orthogonal basis, in principle, it can be more localized than the D-wavelets. While it is true that the biorthogonal wpvb basis is delocalized, by exchanging the roles of the wpvN and wpvb bases and using the wpvN functions to calculate the coefficients, potentially higher compression factors can be achieved than with D-wavelets. Furthermore, there is no need for “mother” and “father” wavelets as in D-wavelets, since there is no need to orthogonalize the high and low resolution functions; there is just a single Gaussian parent wavelet. (3) The wpvN and wpvb bases at critical sampling span exactly the same Hilbert space as the FG method. This exact equivalence with the FG method goes against the conventional wisdom that the sampling density of a windowed Fourier transform has to be larger than the Nyquist density in order to ensure completeness of the basis set (or frame).22 In our experience, no stability problems have been encountered, either before or after discarding basis functions.

In summary, we have presented an extension of the pvb formalism, where the centers and widths of the Gaussian functions obey a scaling relation. This scaled version, which we call wpvb, allows the basis functions to match the classical phase space contours. The wpvb method is particularly useful when the classical phase space area is not “smooth” in the sense that there are sharp peaks in p, as in Coulomb systems. As shown in Ref. 14, the savings are expected to grow dramatically with dimensionality. We therefore believe that the wpvb method may provide an efficient alternative to existing methods for performing electronic structure calculations on atoms and molecules, particularly for excited states. Tests along these lines are in progress.

1.
S. G.
Mallat
,
IEEE Trans. Pattern Anal. Mach. Intell.
11
,
674
(
1989
).
2.
I.
Daubechies
,
IEEE Trans. Inf. Theory
36
,
961
(
1990
).
3.
G. K.
Wallace
,
Commun. ACM
34
,
30
(
1991
).
5.
D.
Sinha
and
A. H.
Tewfik
,
IEEE Trans. Signal Proc.
41
,
3463
(
1993
).
6.
S. P.
Srinivason
and
L. H.
Jamieson
,
IEEE Trans. Signal Proc.
46
,
1085
(
1998
).
7.
D.
Permann
and
I.
Hamilton
,
J. Chem. Phys.
100
,
379
(
1994
).
8.
K.
Cho
,
T. A.
Arias
, and
J. D.
Joannopoulos
,
Phys. Rev. Lett.
71
,
1808
(
1993
).
9.
J. P.
Modisette
,
P.
Nordlander
,
J. L.
Kinsey
, and
B. R.
Johnson
,
Chem. Phys. Lett.
250
,
485
(
1996
).
10.
B.
Poirier
and
A.
Salam
,
J. Chem. Phys.
121
,
1690
(
2004
).
11.
A. I.
Neelor
and
S.
Goedecker
,
J. Comput. Phys.
217
,
312
(
2006
)
12.
J.
von Neumann
,
Math. Ann.
104
,
570
(
1931
).
13.
D.
Gabor
,
J. Inst. Electr. Eng.
93
,
429
(
1946
).
14.
A.
Shimshovitz
and
D. J.
Tannor
,
Phys. Rev. Lett.
109
,
070402
(
2012
).
[PubMed]
15.
R.
Kosloff
in
Numerical Grid Methods and their Application to Schrödinger's Equation
, edited by
C.
Cerjan
(
Kluwer
,
Boston
,
1993
).
16.
C. C.
Marston
and
G. G.
Balint-Kurti
,
J. Chem. Phys.
6
,
3571
(
1989
).
17.
D. T.
Colbert
and
William H.
Miller
,
J. Chem. Phys.
96
,
1982
(
1992
).
18.
D. J.
Tannor
,
Introduction to Quantum Mechanics: A Time-Dependent Perspective
(
University Science Books
,
Sausalito
,
2007
), eq.11.163.
19.
S. R.
Dooley
and
A. K.
Nandi
,
IEEE Trans. Signal Proc.
48
,
1201
(
2000
).
20.
Reference 18, eq.11.172.
21.
J.
Morlet
,
G.
Arens
,
I.
Fourgeau
, and
D.
Giard
,
Geophys.
47
,
203
(
1982
).
22.
I.
Daubechies
,
Ten Lectures on Wavelets
,
CBMS-NSF Regional Conference Series in Applied Mathematics Vol. 61
(
SIAM
,
Philadelphia
,
1992
).