1. Cylindrical DVR Basis¶

Here we discuss the DVR discretization for systems with cylindrical symmetry. For this, one can use the $d=2$ Bessel function basis. Here we specifically discuss and demonstrate how to work with this particular basis as it demonstrates a few potential gotchas. First we give a very brief review of the DVR basis and then we demonstrate the cylindrical basis. $\newcommand{\d}{\mathrm{d}}\newcommand{\vect}[1]{\vec{#1}}\newcommand{\abs}[1]{\lvert#1\rvert}\newcommand{\I}{\mathrm{i}}\newcommand{\braket}[1]{\langle#1\rangle}\newcommand{\ket}[1]{\left|#1\right\rangle}\newcommand{\bra}[1]{\left\langle#1\right|}\newcommand{\pdiff}[2]{\frac{\partial #1}{\partial #2}}\newcommand{\diff}[3][]{\frac{\d^{#1}#2}{\d{#3}{}^{#1}}}$

Table of Contents¶

• 1. Cylindrical DVR Basis
• 2. DVR Basis
  • 2.1 Bessel Function Basis
  • 2.2 Harmonic Oscillator
• 3. Laplacian
• 4. Structure of SLDA code
• 5. Evaluating Bessek Functions

2. DVR Basis¶

The idea of a discrete variable representation (DVR) basis is to find a consistent set of orthonormal basis functions $F_n(r)$ and abscissa $r_m$ such that is quasi-local:

$$ F_n(r_m) = \frac{\delta_{mn}}{\sqrt{w_n}}, \qquad f(r) = \sum_n f_nF_n(r), \qquad f_n = \sqrt{w_n}f(r_n). $$

The quasi-local property of the basis functions means that expanding an arbitrary function in the basis is a simple as evaluating it at the abscissa as shown above. The interpretation of the numbers $w_n$ is that they act as integration quadrature weights which are exact for the product $g^*(r)f(r)$ of any functions represented exactly in the basis:

$$ \int g^*(r) f(r) = \sum_{mn} \sqrt{w_mw_n}g^*(r_m)f(r_n) \int F_m(r)^*F_n(r) = \sum_{n} w_n g^*(r_n)f(r_n). $$

These basis functions are typically $\delta$ functions projected onto a space of limited momenta with $k<k_c$ and are often constructed to deal with particular boundary conditions. In particular, here we consider the Bessel function basis which is constructed to represent the radial wavefunction of centrally symmetric potentials.

One can often express the kinetic energy of a Hamiltonian in this basis as a matrix $K_{mn} = \braket{F_m|\op{K}|F_n}$ and then obtain exponential accuracy with a Hamiltonian whose potential is diagonal $V_{mn} = \delta_{mn} V(r_n)$ and obtained by simply evaluating the potential at the abscissa. This only works if the potential is smooth (otherwise high momenta $k > k_c$ appear). Thus, to deal with the singularity in the effective central potential, one needs to construct the basis explicitly to incorporate this into the "kinetic" energy. This is the purpose of the Bessel function basis described below.

2.1 Bessel Function Basis¶

For central problems one expands the radial wavefunction $u(r) = r^{(d-1)/2}\psi(r)$ for each angular momentum $\nu = l + d/2 - 1$. The Bessel function DVR basis described here provides a basis for smooth radial wavefunctions, explicitly treating the singular centrifugal piece as part of the kinetic energy:

$$ K_{mn} = \braket{F_m|\op{K}_{\nu}|F_n} = \braket{F_m|\left(-\I\pdiff{}{r}\right)^2 + \frac{\nu^2 - \tfrac{1}{4}}{r^2}|F_n}, \\ u(r) = \sum_{n} \sqrt{w_n}u(r_n)F_n(r) = \sum_{n} \sqrt{w_n}r_n^{(d-1)/2}\psi(r_n)F_n(r),\\ \psi(r) = \frac{u(r)}{r^{(d-1)/2}} = \sum_{n} \sqrt{w_n} r_n^{(d-1)/2}\psi(r_n)\frac{F_n(r)}{r^{(d-1)/2}}. $$

In principle, there is a different basis for each angular quantum number $l$, but in practice, one can often use just a few basis: in $d=3$ for example, one can use the $l=0$ basis for all even $l$ and the $l=1$ basis for all odd $l$. The abscissa and basis functions are defined in terms of the Bessel functions $J_\nu(z)$ and their roots $z_{\nu n}$ (see~\cite{LC:2002}):

$$ \begin{gather} \diff[2]{u(r)}{r} - \frac{\nu^2 - 1/4}{r^2}u(r) + \frac{2m}{\hbar^2}[E - V(r)]u(r) = 0\\ \begin{aligned} J_\nu(z_{\nu n}) &= 0, & w_n &= \frac{2}{k_c z_{\nu n}J'_{\nu}(z_{\nu n})^2}, & F_n(r) &= (-1)^{n+1}\frac{k_c z_{\nu n}\sqrt{2r}}{k_c^2r^2 - z_{\nu n}^2} J_\nu(k_c r) \end{aligned}\\ \begin{aligned} K_{m\neq n} &= \frac{8k_c^2(-1)^{m-n}z_{\nu n}z_{\nu m}}{(z_{\nu n}^2 - z_{\nu m}^2)^2}, & K_{nn} &= \frac{k_c^2}{3}\left[1 + \frac{2(\nu^2 - 1)}{z_{\nu n}^2}\right] \end{aligned} \end{gather} $$

2.2 Harmonic Oscillator¶

We start by solving the 2D Harmonic Oscillator.

Here is a simple class that will define the DVR basis in terms of a number of abscissa N and momentum scale k.

Here we demonstrate how well the basis works. For the explicit l values, the errors are very small. When we fix l_explicit = 0, then even l work well, but the odd l are not so good. Likewise, when we fix l_explicit = 1, then the odd l work well, but not the even l. Thus, we have a similar situation to the 3D case where we should use the abscissa for l=0 and l=1, and then represent the higher momenta by adding the appropriate correction as a potential.

3. Laplacian¶

Another application of the DVR basis is to compute the Laplacian by acting the kinetic energy operator $-\op{K}$ on states. Here we demonstrate this with a Gaussian in $d$ dimensions:

$$ \nabla^2 e^{-(r/r_0)^2/2} = \frac{r^2 - dr_0^2}{r_0^4}e^{-(r/r_0)^2/2} $$

For the cylindrical basis, $d=2$.

Here we must be a little careful to recall that the basis is a basis for the radial wavefunctions and one must not forget the weights. Thus:

$$ -\nabla^2 \Psi(\vect{r}) = -\nabla^2 Y_{lm}(\Omega)\psi_l(r) = \frac{Y_{lm}(\Omega)}{r^{(d-1)/2}} \op{K}_{l} u_l(r). $$

In particular, for fully summetric functions $\psi(r)$ (i.e. $l=0$), the Laplacian is computed as:

$$ -\nabla^2 \psi(r) = \sum_{mn} \frac{F_{m}(r)}{\sqrt{w_m}r^{(d-1)/2}} K_{mn} \left(\sqrt{w_n}r_n^{(d-1)/2}\psi(r_n)\right) $$

4. Structure of SLDA code¶

The general structure of the time-dependent SLDA code is as follows:

• Define initial state.

• Distribute wavefunctions (u_n and v_n) amoungst the nodes

5. Evaluating Bessel Functions¶

This section contains some notes about numerically evaluating the Bessel functions, especially near the roots when there is a pole in the denominator:

$$ \frac{\sqrt{z}J_{\nu}(z)}{z - z_{n}} $$

As $z$ approaches $z_n$, this has the form of 0/0, so one can apply a form of l'Hopital's rule to reduce the round-off error. The specified form of the function has been chosen for special properties of the Bessel functions. Express the function as

\begin{align} F(z) &= \frac{f(z)}{z - z_n} = \frac{\sqrt{z}J_{\nu}(z)}{z - z_n}\\ F'(z) &= \frac{f'(z)}{z - z_n} - \frac{f(z)}{(z - z_n)^2} \end{align}

Let $\delta = z - z_n$. Close to the singular point we use the Taylor series:

$$ \sum_{m=0}^{\infty}\frac{a_m\delta^{m}}{m!} $$\begin{align} H(z) &= h'(z_n) + \sum_{m=3}^{\infty}\frac{h^{(m)}(z_n)\delta^{m-1}}{m!} = \sum_{m=0}^{\infty}\frac{h^{(m+1)}(z_n)\delta^{m}}{(m+1)m!}\\ a_m &= \frac{h^{(m+1)}}{m+1}\\ H'(z) &= \sum_{m=3}^{\infty}\frac{(m-1)h^{(m)}(z_n)\delta^{m-2}}{m!} = \sum_{m=1}^{\infty}\frac{h^{(m+2)}(z_n)\delta^{m}}{(m+2)m!}\\ a_m &= \frac{h^{(m+2)}}{m+2} \end{align}

The first few derivatives are presented here:

\begin{align} h(z) &= \sqrt{z}J_{\nu}(z)\\ h'(z) &= \frac{J_{\nu}(z)}{2\sqrt{z}} + \sqrt{z}J'_{\nu}(z) = \frac{f(z)}{2z} + \sqrt{z}J'_{\nu}(z)\\ h''(z) &= \sqrt{z}J_{\nu}(z)\left( \frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right) = f(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)\\ h'''(z) &= f'(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right) - 2f(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3}\\ h^{(4)}(z) &= h(z)\left[ \left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)^2 + 6\frac{\nu^2 - \tfrac{1}{4}}{z^4}\right] - 4f'(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3} \end{align}

(These were checked with Maple.)

Evaluated at the root $z=z_n$ these become:

\begin{align} h(z_{n}) &= 0\\ h'(z_{n}) &= \sqrt{z_{n}}J'_{\nu}(z_{n})\\ h''(z_{n}) &= 0\\ h'''(z_{n}) &= h'(z_{n})\left( \frac{\nu^2 - \tfrac{1}{4}}{z_{n}^2} - 1\right)\\ h^{(4)}(z_{n}) &= - 4h'(z_{n})\frac{\nu^2 - \tfrac{1}{4}}{z_{n}^3} \end{align}

with both the function and the second derivative vanishing.

To determine where to use this formula, we match the estimate roundoff error with the truncation error. The Bessel functions are of order unity and are typically calculated to an absolute accuracy of $\epsilon$. The round-off error in the numerator is $\epsilon f(z)$ and $\epsilon \sqrt{2}z_n$ in the denominator. The roundoff errors in the denominator dominate both cases:

\begin{align} \delta H(z) &\sim \epsilon \frac{\sqrt{2}z_n H(z)}{\delta} \sim \frac{\sqrt{2}\epsilon z_n h(z)}{\delta^2} \sim \frac{\sqrt{2}\epsilon z_n h'(z_n)}{\delta}\\ \delta H'(z) &\sim \frac{2\epsilon z_n h(z)}{\delta^3} \sim \frac{2\epsilon z_n h'(z_n)}{\delta^2} \end{align}

To choose the appropriate transition point, we equate half of this with the truncation error to transition points:

\begin{align} \delta_c &\sim \left( \frac{72\epsilon z_n h'(z_n)}{\sqrt{2}h^{(4)}(z_n)} \right)^{1/4} \sim \left(\frac{72\epsilon z_n}{\sqrt{2}} \right)^{1/4}\\ \delta_c' &\sim \left(120\epsilon z_n\right)^{1/5} \end{align}

the fact that $h(z)$ behaves asymptotically as a $\sqrt{2/\pi}\cos(z + \phi)$ and so all derivatives have essentially the same magnitude.