Numerics.TISE

Solving the time-independent Schrödinger equation (TISE)

Conventional approaches to solve the time-independent Schrödinger equation

are based on a (exact or approximate) diagonalization of a suitably chosen matrix representation of the Hamiltonian operator. This is implemented in WavePacket using different discrete variable representation (DVR) schemes,. Basically, we expand the wave function that is sought in a set of functions (DVR grid functions), represent the Hamiltonian as a matrix in this set, and diagonalize the resulting matrix, where the (optional) use of sparse matrix techniques is crucial for multidimensional applications. When choosing an FFT based DVR scheme, this approach is equivalent to the popular Fourier Grid Hamiltonian (FGH) method by Meyer, later rediscovered by Balint-Kurti..

As an alternative to direct diagonalization, propagation in imaginary time relaxes an initial wavefunction toward the ground state.

Direct diagonalization

The most direct method for solving the time-independent Schrödinger equation is the setup of a DVR-based matrix representation of the underlying Hamilton operator, and the diagonalization of the resulting matrix. This gives accurate results even for highly excited states, but the matrices tend to get very large for high-dimensional problems (but sparsity can often be exploited). Also, for almost degenerate eigenstates, symmetries may not be preserved due to numerical errors. The diagonalization scheme is available for all potentials, employing the "natural" kinetic energy operators of the various grids (FGH method or other DVRs based on angular momentum algebra).

One of the best-known methods here is the Fourier Grid Hamiltonian (FGH) method for plane wave expansions. It was originally discovered by R.Meyer, later rediscovered by C.C.Marston and G.G.Balint-Kurti, and also generalized to the multidimensional case (see Eq. (13) in work by S. P. Webb and Sh. Hammes-Schiffer)

A comparison of the FGH method with a Gauss-Hermite expansion can be found in the Double well demo. A direct diagonalization in a spherical harmonics basis is given in the Pendular states demo

Imaginary time propagation (relaxation method)

Wavepacket propagation in "imaginary time" t→-iτ transforms the TISE into the form of a diffusion equation with time evolution operator

where the normalized Hamiltonian (spectrum between -1 and +1) is defined as

where E_min and ΔE are the minimum energy and the spectral range of the Hamiltonian, respectively, and where the dimensionless "Kosloff parameter" is defined as

It is then advantageous to expand the imaginary time evolution operator in a series of real Chebychev polynomials

where the coefficients c_n are modified Bessel functions (which decay quickly) and where the φ_n are the real Chebychev polynomials which are calculated conveniently using the recursion

The advantage is that if n large enough, the Bessel functions decay fast to zero! In practice, it has turned out that α shouldn't be below 10.

Since the norm of the ground states decays slower than that of the excited states, this propagation converges towards the ground state of a quantum system. However, the Chebychev method can also be used to obtain (at least the lowest few) excited states. In order to calculate the N-th bound state, we apply the operator (1-P) H (1-P) instead of H in the Chebychev scheme, where P is the projection on the space spanned by the lower N-1 states. This procedure can be repeated projecting out any ground state contribution to yield the first excited state, etc. In practice, the propagation time required for convergence depends on the energy gap between the ground and first excited states. Also, when trying to find very highly excited states, errors quickly grow out of hand, so only a handful of excited states can be reliably calculated.

• see original work by R. Kosloff and H. Tal-Ezer
• See the Morse oscillator demo for a demonstration with one of our toy systems
• See the H3+ molecule demo for a demonstration with a real molecule