# Non-hermitian discretizations in quantum mechanics

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?

• I don't know anything about quantum mechanics or solving Schroedinger equation, but given 2nd derivative matrix should theoretically be symmetric, wouldn't you usually symmetrize the approximated 2nd derivative matrix prior to use? That's what I do when get slight asymmetries in calculated (automatic differentiation or finite differences) 2nd derv. matrix due to roundoff. If H is the nonsymmetric 2nd derivative matrix, then use $0.5*(H + H^T)$ or replace 1 triangle by transpose of the other. Or just use the real part of the eigenvalues. If your asymmetry is extreme, is anythiig meaningful? – Mark L. Stone May 24 '18 at 23:01
• Substitute Hermitian and conjugate transpose as appropriate in my above comment. – Mark L. Stone May 24 '18 at 23:44
• That the matrix isn't Hermitian does not imply that the eigenvalue are not real. – nicoguaro May 25 '18 at 4:26
• @MarkL.Stone: thanks for your comment. I know this solution, but I always have the feeling it spoils the effect of the non-equidistant grid. Consider two gridpoints $x_i$ and $x_{i+1}$ and a three point stencil. Symmetry means, the contribution of point $x_{i}$ to the derivative at point $x_{i+1}$ is the same as the contribution to contribution of point $x_{i+1}$ to the derivative at point $x_{i}$. But, depending on the other points $x_{i-1}$ and $x_{i+2}$, the one can be significant and the other negligible. – davidhigh May 25 '18 at 16:49

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

• Thanks for your answer. Can you give some more details? I see that there can be examples where one has only real eigenvalues. In general, however, one has those complex eigenvalues, and those completely spoil the usual concepts (e.g. the time propagation by exponentiation of the eigenvalues). When exactly have non-hermitian matrices which discretize self-adjoint operators only real eigenvalues? – davidhigh May 25 '18 at 23:45
• @davidhigh, I don't have a theoretical answer off the top of my head. Nevertheless, every time I used non-Hermitian discretizations for problems with real eigenvalues (mostly frequencies and energies) I obtained real eigenvalues (maybe with negligible imaginary parts). I have obtained these results with finite differences, collocation methods and coupled finite elements. – nicoguaro May 25 '18 at 23:55