Non-hermitian discretizations in quantum mechanics

Question

Consider the Schroedinger equation

$$\left(-\frac12\frac{\partial^2}{\partial x^2} + V(x) \right) \psi(x) = E \psi(x)$$

The usual way to solve it is to introduce a discretization of $\psi(x)$. This yields an eigenvalue problem which can be solved by standard methods from linear algebra. When using a discretization that leads to hermitian (resp. symmetric) matrices, one obtains real eigenvalues which determine the energies of the eigenstates of the system.

However, when one uses a non-symmetric discretization for the second-derivative matrix, such as those from a Chebyshev or Legendre collocation, or simply from an equidistant grid with forward or backward differences at the boundaries, one in general obtains complex eigenvalues.

Does that mean one can't use non-symmetric discretizations for quantum mechanics? Are there any workarounds?

Answer 1

If your discrete equations represent your problem properly you should have the right eigenvalues, although numerically you can get a small imaginary part when solving the eigenvalue problem.

As an example, let us consider the axisymmetric case of a cylindrical well. The Schrödinger equation would be written as

$$\frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r}\frac{\partial \psi}{\partial r} = E \psi\, .$$

Using a naive finite difference we can use the following approximations

\begin{align} &\frac{\partial^2 \psi(r_i)}{\partial r^2} \approx \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{\Delta r^2}\, ,\\ &\frac{\partial \psi(r_i)}{\partial r} \approx \frac{\psi_{i+1} - \psi_{i-1}}{2\Delta}\, , \end{align}

where the second approximation leads to a skew-symmetric matrix. The discrete problem would be then

$$[D_{2r} + R_\text{inv} D_{r}] \{\Psi\} = E \{\Psi\}\, ,$$

being $D_{2r}$ the matrix of second derivatives, $D_{r}$ the matrix of first derivatives and $R_\text{inv}$ a diagonal matrix with $\frac{1}{r_i + \Delta r/2}$ as entries.

The following snippet solves this problem

from __future__ import division, print_function
from scipy.special import j0, jn_zeros
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigs
import matplotlib.pyplot as plt

# Parameters
r_max = 1.0
dr = 0.01
r = np.arange(0, r_max, dr)
npts = r.shape[0]
nvals = 10

# Matrices
D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts - 1, npts - 1))/dr**2
D1 = diags([-0.5, 0, 0.5], [-1, 0, 1], shape=(npts - 1, npts - 1))/dr
r_new = r[0:npts-1] + dr/2
R_inv = diags(1/r_new, 0)
A = D2 + R_inv.dot(D1)

# Finite differences
vals, vecs = eigs(-A, k=nvals, which="SM") 
vecs = np.vstack((vecs, np.zeros(vecs.shape[1])))

# Analytic results
vals_anal = jn_zeros(0, nvals)
vecs_anal = j0(np.outer(r, vals_anal)/r_max)


## Plots

# Eigenvalues
plt.figure()
plt.plot(vals_anal**2)
plt.plot(vals, "ok")
plt.legend(["Analytic eigenvalues",
            "Numeric eigenvalues"])
plt.xlabel(r"$n$")
plt.ylabel(r"$E_n$")


# Eigenvectors
plt.figure()
plt.subplot(121)
plt.plot(r, np.abs(vecs[:, 0:4]/vecs[0, 0:4])**2)
plt.title("Finite differences")
plt.xlabel(r"$r$")
plt.ylabel(r"$|\psi|^2$")

plt.subplot(122)
plt.plot(r, np.abs(vecs_anal[:, 0:4])**2)
plt.title("Analytic solution")
plt.xlabel(r"$r$")
plt.ylabel(r"$|\psi|^2$")

plt.tight_layout()
plt.show()

With the following eigenvalues

and these eigenvectors