Spectral Methods: Singular Potentials

Spectral Methods: Singular Potentials#

To obtain this order of accuracy, we must be free to choose both the weights and the points. If we want \(x_0\) and/or \(x_{N-1}\) as the endpoints of the integration interval, then we must sacrifice one or two orders:

Gauss-Radau rules or Radau quadrature fixes one endpoint and is exact up to order \(2N-2\).
Gauss-Lobatto rules or Lobatto quadrature fixes both endpoints and is exact up to order \(2N-3\).

On the surface, this sees to be a problem for constructing the DVR basis, but as pointed out in [Schneider and Nygaard, 2004], the Radau quadrature based on the Laguerre polynomials is find for solving the radial equation if we exclude the first abscissa \(x_0 = 0\) where \(u(0) = 0\). Thus, we can use an \(N+1\)-point Radau quadrature, exact for polynomials of order \(2N\), and take the \(N\) positive abscissa and weights to construct the DVR basis.

Laguerre-Radau Quadrature#

The generalized Laguerre polynomials \(L_{N}^{(\alpha)}\) are orthogonal on \([0, \infty)\) with weight \(w(x) = x^{\alpha}e^{-x}\). The Radau quadrature abscissa \(x_i\) include \(x_0=0\) and the roots of \(L_{N}^{(\alpha+1)}(x)\), with quadrature weights \(w_i\) [Gautschi, 2000]:

\[\begin{gather*} L_{N}^{(\alpha+1)}(x_i) = 0, \qquad w_{0} = \frac{\Gamma(\alpha+1)}{{N+\alpha+1} \choose {N}}, \qquad w_{i} = \frac{\Gamma(\alpha+1)}{N+\alpha+1} {N+\alpha \choose N} \frac{1}{[L_{N}^{(\alpha)}(x_i)]^2}. \end{gather*}\]

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 8
      4 
      5 import numpy as np
      6 from scipy.integrate import quad
      7 
----> 8 from phys_581.dvr import roots_genlaguerre_radau, LaguerreRadauDVR
      9 
     10 dvr = LaguerreRadauDVR(N=10)
     11 ns = range(dvr.N+1)

ModuleNotFoundError: No module named 'phys_581'

Kinetic Energy#

To compute the kinetic energy we note that the basis functions are a sum of terms of the form

\[\begin{gather*} f(r) = \overbrace{\sqrt{r^{\alpha}e^{-r}}}^{\sqrt{w}}L(r),\qquad f'(r) = \sqrt{w(r)}\left(\frac{\alpha - r}{2r}L + L'(r)\right). \end{gather*}\]

Thus, we can compute the kinetic energy matrix using the quadrature forumlae:

\[\begin{gather*} K_{mn} = \braket{u_m|\left(-\tfrac{1}{2}\diff[2]{}{r}\right)|u_n} = \frac{1}{2}\int\d{r}\; u_{m}'(r) u_{n}'(r)\d{r} = \frac{1}{2}\sum_{i} w_i \frac{u_{m}'(r_i) u_{n}'(r_i)}{w(r_i)} \end{gather*}\]

Detailed expressions for \(N=1\) (used to check code).

Getting everything working here is quite challenging as there are lots of places for making mistakes. Ultimately we should have a comprehensive set of tests, but initially it is helpful to work with a small basis. Working through these and explicitly checking the code helped me find several subtle issues. We take \(N=1\) with two abscissa.

With \(\alpha = 0\), we have explicitly:

\[\begin{align*} P_0(r) = L_0^{(\alpha)}(r) &= 1\\ P_1(r) = L_1^{(\alpha)}(r) &= 1 +\alpha - r = 1-r\\ L_{1}^{(\alpha+1)}(r) &= 2+\alpha-r = 2 - r\\ w_0 = w_1 &= \frac{1}{2},\\ \mat{C} &= \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1\\ 1 & -1 \end{pmatrix},\\ F_0(r) &= \frac{P_0(r) + P_1(r)}{\sqrt{2}} = \frac{2-r}{\sqrt{2}}\\ F_1(r) &= \frac{P_0(r) - P_1(r)}{\sqrt{2}} = \frac{-r}{\sqrt{2}}\\ u_1(r) &= \tfrac{-1}{\sqrt{2}}re^{-r/2}. \end{align*}\]

The abscissa are \(r_0 = 0\) and \(r_1 = 2\), the latter being the roots of \(L_{N}^{(\alpha+1)}(r)\).

For \(\alpha=1\) we have

\[\begin{align*} P_0(r) = L_0^{(1)}(r) &= 1\\ P_1(r) = \tfrac{1}{\sqrt{2}}L_1^{(1)}(r) &= \frac{2 - r}{\sqrt{2}}\\ L_{1}^{(2)}(r) &= 3 - r\\ (w_0, w_1) &= (\tfrac{1}{3}, \tfrac{2}{3}),\\ \mat{C} &= \frac{1}{\sqrt{3}} \begin{pmatrix} 1 & \sqrt{2}\\ \sqrt{2} & -1 \end{pmatrix},\\ F_0(r) = \frac{P_0(r) + \sqrt{2}P_1(r)}{\sqrt{3}} &= \frac{3 - r}{\sqrt{3}},\\ F_1(r) = \frac{\sqrt{2}P_0(r) - P_1(r)}{\sqrt{3}} &= \frac{r}{\sqrt{6}},\\ u_0(r) &= \tfrac{1}{\sqrt{3}}(3-r)\sqrt{r}e^{-r/2},\\ u_1(r) &= \tfrac{1}{\sqrt{6}}r\sqrt{r}e^{-r/2} \end{align*}\]

Thus, despite the fact that the \(\alpha = 1\) basis has \(L_{n}^{(1)}\) basis polynomials, the factor of \(\sqrt{r}\) ruins the convergence.

Note that for \(m = \alpha = \hbar = 1\), the first few state for Hydrogen are

\[\begin{gather*} u_{01}(r) = 2re^{-r}L_{0}^{(1)}(2r) = 2re^{-r}, \\ u_{02}(r) = \tfrac{1}{2\sqrt{2}} re^{-r/2}L_{1}^{(1)}(r) = \tfrac{1}{4}e^{-r/2}(2-r)r. \end{gather*}\]

One final comment: We have expressed our basis in some sort of “natural” units where \(r\) is dimensionless. This is no good for physics. Instead, we must choose a length scale \(a\) in which to express our problem. Once we choose, this scale, the basis can be used by scaling \(r\rightarrow r/a\) to be dimensionless, and then adding the missing factors \(\mat{K}/a^2\) to the kinetic energy. To choose \(a\), note that the quadrature is exact for radial wavefunctions of the form (setting \(\alpha = 0\) here)

\[\begin{gather*} u(r) \propto e^{-r/2a}L\left(\frac{r}{a}\right). \end{gather*}\]

This might indicate that an optimal choice is \(a \propto n \propto 1/\sqrt{-E}\) for Coulomb-like potentials. We will explore how well this works below, but note that choosing a slightly incorrect \(a\) still works well because

\[\begin{gather*} e^{-r/2a + \epsilon r} = e^{-r/2a}\left(1 + \epsilon r + \frac{\epsilon^2}{2!}r^2 + \cdots\right). \end{gather*}\]

Thus, the errors induce polynomial factors that are well represented in the basis.

Example: Coulomb-like Potentials.#

Here we use the Laguerre-Radau approach described above to solve for states in a Coulomb-like potential.

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 15
     11 ns = 1+np.arange(Nstates)
     12 E_exact = -1/2/ns**2
     13 
     14 Ns = np.arange(Nstates, Nmax+1, 4)
---> 15 errs = [get_Es(N, a=a)[:Nstates]/E_exact - 1 for N in Ns]
     16 
     17 fig, ax = plt.subplots()
     18 ax.semilogy(Ns, np.abs(errs), '-+')

Cell In[3], line 2, in get_Es(N, a)
      1 def get_Es(N=32, a=1.0):
----> 2     dvr = LaguerreRadauDVR(N=N, alpha=0)
      3     H = dvr.T[1:, 1:] + np.diag(-a/dvr.r[1:])
      4     En = np.linalg.eigvalsh(H)/a**2
      5     return En

NameError: name 'LaguerreRadauDVR' is not defined

Here we see the convergence of the lowest 15 states as a function of basis size \(N\). Around \(N \approx 180\), numerical errors in the computation start to corrupt the basis, leading to zero weights (see warnings), overflows, etc. While there might be ways to mitigate this (we have included one in the code, using exact integer math for some of the coefficients), it is probably not worth pursuing this direction.

Instead, it us much more profitable to explore adjusting the length scale \(a\).

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 8
      4 ns = 1+np.arange(Nstates)
      5 E_exact = -1/2/ns**2
      6 
      7 Ns = np.arange(Nstates, Nmax+1, 4)
----> 8 errs = [get_Es(N, a=a)[:Nstates]/E_exact - 1 for N in Ns]
      9 
     10 fig, ax = plt.subplots()
     11 ax.semilogy(Ns, np.abs(errs), '-+')

Cell In[3], line 2, in get_Es(N, a)
      1 def get_Es(N=32, a=1.0):
----> 2     dvr = LaguerreRadauDVR(N=N, alpha=0)
      3     H = dvr.T[1:, 1:] + np.diag(-a/dvr.r[1:])
      4     En = np.linalg.eigvalsh(H)/a**2
      5     return En

NameError: name 'LaguerreRadauDVR' is not defined

Here we will present without proof a spectral method based on the Bessel-function discrete variable representation (DVR). Here we introduce a basis \(\ket{F_{\nu,n}}\) obtained by projecting onto a space with wave-vectors less than \(\abs{k} < k_\max\):

\[\begin{gather*} \op{P} = \int_{\abs{\vect{k}}<k_\max}\!\!\!\!\!\d^{d} \vect{k}\; \ket{\vect{k}}\bra{\vect{k}}. \end{gather*}\]

In this basis, the following representation for the operator \(\op{K}\) is exact:

\[\begin{gather*} \mat{K}^{(\lambda)}_{m,n} = \braket{F_{\nu, m}|\left(-\diff[2]{}{r} + \frac{\nu^2 - 1/4}{r^2}\right)|F_{\nu, n}}\\ = k_{\max} \begin{cases} \frac{1}{3}\left(1 + \frac{2(\nu^2 - 1)}{z_{\nu,n}^2}\right) & m=n,\\ (-1)^{n-m}\frac{8 z_{\nu,m}z_{\nu,n}}{(z_{\nu, m}^2 - z_{\nu, n}^2)^2}, & m \neq n, \end{cases}\\ J_{\nu}(z_{\nu, n}) = 0, \end{gather*}\]

where \(z_{\nu, n}\) are the roots of the Bessel functions of the first kind.

Furthermore, the basis is quasi-local, so that the potential operator can be expressed as a diagonal matrix

\[\begin{gather*} \braket{F_{\nu, m}|V(\op{r})|F_{\nu, n}} \approx V(r_{\nu, n}) \delta_{mn}, \qquad r_{\nu, n} = \frac{z_{\nu, n}}{k_\max}. \end{gather*}\]

This is not exact, but provides exponential accuracy for analytic potentials.

Basis functions

We will not need them for our work here, but for reference, the basis functions are:

\[\begin{gather*} \braket{F_{\nu, n}|F_{\nu, m}} = \int_0^{\infty}\d{r}\; F_{\nu, n}(r)F_{\nu, m}(r) = \delta_{mn},\\ F_{\nu, n}(r) = (-1)^{n+1} \frac{r_{\nu, n} \sqrt{2r}} {r^2 - r_{\nu, n}^2}J_\nu(k_\max r). \end{gather*}\]

Note

One must be careful evaluating these near the roots where the denominator vanishes. In our code, we do this with a careful application of Taylor series and L’Hopital’s rule.

Note that these are zero at all abscissa \(r_{\nu, m}\) except for \(m = n\):

\[\begin{gather*} F_{\nu, n}(r_{\nu, m}) = \delta_{n,m} \underbrace{ \sqrt{k_\max} (-1)^{n+1} \sqrt{\frac{z_{\nu, n}}{2}} J_\nu'(z_{\nu, n}) }_{w_n}. \end{gather*}\]

This is a key property of DVR bases which are closely related to the classical orthogonal polynomials where a careful choice of both abscissa and weights leads to twice the expected accuracy when integrating. It also allows us to express a function in the basis by simply computing its value at the abscissa:

\[\begin{gather*} f(r) = \sum_{n} f_n F_{\nu, n}(r), \qquad f(r_{\nu, m}) = \sum_{n} f_n F_{\nu, n}(r_{\nu, m}) = f_m w_{m},\\ f_n = \braket{F_{\nu, n}|f} = \frac{f(r_{\nu, n})}{w_n}. \end{gather*}\]

Finally, we note that from these properties, the coefficients \(w_n\) act as integration weights for functions expressed in the basis:

\[\begin{gather*} \int_0^{r} f(r) \d{r} = \sum_{n} f_{n} \int_0^{r} F_{\nu, n}(r) \d{r} = \end{gather*}\]

Demonstration#

As a quick demonstration, we find the eigenstates of the spherical harmonic oscillator, and the eigenstates of hydrogenic atoms:

\[\begin{gather*} V_{HO}(r) = \frac{m\omega^2r^2}{2}, \qquad E_{l,n} = \hbar\omega\left(2n + l + \frac{d}{2}\right), \qquad a_{ho} = \frac{\hbar}{\sqrt{m\omega}},\\ V_{H}(r) = \frac{-\alpha}{r}, \qquad E_{l,n} = \underbrace{ -\frac{m \alpha^2}{2\hbar^2}}_{-13.6\;\mathrm{eV}} \frac{1}{(l+n+1)^2},\qquad a_{h} = \frac{\hbar^2 }{m \alpha}. \end{gather*}\]

The numerical value is given for hydrogen where \(\alpha = e^2/4\pi\epsilon_0 \approx 14.4\)eV Å.

from phys_581 import bessel, dvr

N = 10
d = 3
hbar = m = w = 1.0
a_ho = hbar / np.sqrt(m*w)
R = N*a_ho
k_max = N/a_ho

def V(r):
    return m * (w*r)**2 / 2
    
def get_E(l, N=N):
    n = np.arange(N)
    return hbar * w * (2*n + l + d/2)

for d in [2, 3, 4]:
  basis = dvr.SphericalDVRBasis(R=R, d=d, k_max=k_max)
  for l in [0, 1, 2, 3, 4]:
    r = basis.get_rn(l=l)
    H = hbar**2 / 2/ m * basis.get_K(l) + np.diag(V(r))
    assert np.allclose(np.linalg.eigvalsh(H)[:N], get_E(l=l))

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[5], line 1
----> 1 from phys_581 import bessel, dvr
      2 
      3 N = 10
      4 d = 3

ModuleNotFoundError: No module named 'phys_581'

from phys_581 import bessel, dvr

N = 5
d = 3
hbar = m = e = alpha = 1.0

a_h = hbar**2 / m / alpha
R = 10*N*a_h
k_max = 10*N/a_h

def V(r):
    return -alpha / r
    
def get_E(l, N=N):
    n = np.arange(N)
    return -m * alpha**2 / 2 / hbar**2 / (1 + n + l)**2

basis = dvr.SphericalDVRBasis(R=R, d=d, k_max=k_max)
for l in [0, 1, 2, 3, 4]:
    r = basis.get_rn(l=l)
    H = hbar**2 / 2/ m * basis.get_K(l) + np.diag(V(r))
    print(d, l)
    print(np.linalg.eigvalsh(H)[:N])
    print(get_E(l=l))
    #assert np.allclose(np.linalg.eigvalsh(H)[:N], get_E(l=l))

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[6], line 1
----> 1 from phys_581 import bessel, dvr
      2 
      3 N = 5
      4 d = 3

ModuleNotFoundError: No module named 'phys_581'

k = 2*np.pi * np.fft.fftfreq(N, dx)
K = (hbar*k)**2/2/m
A = np.fft.ifft(K*np.fft.fft(np.eye(N), axis=1), axis=1)
V = np.diag(Vx)
C = None

E, psi = sp.linalg.eig(A+V, C)
inds = np.argsort(abs(E))
E, psi = E[inds], psi[:, inds]
E[:10].real

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 k = 2*np.pi * np.fft.fftfreq(N, dx)
      2 K = (hbar*k)**2/2/m
      3 A = np.fft.ifft(K*np.fft.fft(np.eye(N), axis=1), axis=1)
      4 V = np.diag(Vx)

NameError: name 'N' is not defined

Spectral Methods: Singular Potentials

Contents

Spectral Methods: Singular Potentials#

Laguerre-Radau Quadrature#

Kinetic Energy#

Example: Coulomb-like Potentials.#

Demonstration#

See Also#