4. The time-independent Schrödinger equation

The time-independent Schrödinger equation

\[\left(-\frac{1}{2}\nabla^2 + V(\mathbf{r}) \right) \Psi_k(\mathbf{r}) = E_k \Psi_k(\mathbf{r}), \ \ k=1,2,3,\cdots\]

For a spherical symmetric potential \(V(r)\), the eigenfunctions can be taken as

\[\Psi_{n,l,m}(\mathbf{r}) = r^{-1} u_{n,l}(r) Y_{l,m}(\theta, \phi),\]

and the (radial) TISE becomes

\[-\frac{1}{2}\frac{d^2 u_{n,l}(r)}{d r^2}+\frac{l(l+1)}{2 r^2} u_{n,l}(r) + V(r)u_{n,l} = \epsilon_{n,l} u_{n,l}(r).\]

Gauss-Legendre-Lobatto quadrature is defined on \(x \in [-1,1]\). We need to map the grid/Lobatto points \(x_i \in [-1,1]\) into radial points \(r(x): [-1,1] \rightarrow [0, r_{\text{max}}]\). In that case, \(u_{n,l}(r) = u_{n,l}(r(x))\).

The (rational) mapping

\[r(x) = L \frac{1+x}{1-x+\alpha}, \ \ \alpha = \frac{2L}{r_{\text{max}}} \label{x_to_r},\]

has often been used in earlier work. The parameters \(L\) and \(\alpha\) control the length of the grid and the density of points near \(r=0\).

Another option is to use the linear mapping given by

\[r(x) = \frac{r_{\text{max}}}{2}(x+1).\]

In order to use the Gauss-Legendre-Lobatto pseudospectral method, we must first formulate the radial TISE with respecto to \(x\).

By the chain rule we find that for an arbitrary function \(\psi(r(x))\)

\[\begin{split}\frac{d \psi}{dr} &= \frac{1}{\dot{r}(x)} \frac{d \psi}{dx}, \\ \frac{d^2 \psi}{dr^2} &= \frac{1}{\dot{r}(x)^2} \frac{d^2 \psi}{dx^2} - \frac{\ddot{r}(x)}{\dot{r}(x)^3} \frac{d \psi}{dx},\end{split}\]

and the radial TISE becomes

(1)\[-\frac{1}{2} \left( \frac{1}{\dot{r}(x)^2} \frac{d^2 u_{n,l}}{dx^2} - \frac{\ddot{r}(x)}{\dot{r}(x)^3} \frac{d u_{n,l}}{dx} \right) + V_l(r(x)) u_{n,l}(r(x)) = \epsilon_{n,l} u_{n,l}(r(x)),\]

where we have defined \(V_l(r(x)) \equiv V(r(x)) + \frac{l(l+1)}{2 r(x)^2}\). Note that this is, in general, an unsymmetric eigenvalue problem due to the presence of the term \(\frac{\ddot{r}(x)}{\dot{r}(x)^3}\).

In order to reformulate the above as a symmetric eigenvalue problem, we define the new wavefunction \(\phi_{n,l}(x) = \dot{r}(x)^{1/2} u_{n,l}(r(x))\).

Insertion of this into Eq. Link title yields a symmetric eigenvalue problem

(2)\[\left(-\frac{1}{2} \frac{1}{\dot{r}(x)} \frac{d^2}{dx^2} \frac{1}{\dot{r}(x)} + V_l(r(x))+\tilde{V}(r(x)) \right) \phi_{n,l}(x) = \epsilon_{n,l} \phi_{n,l}(x),\]

where we have defined

\[\tilde{V}(x) \equiv \frac{2\dddot{r}(x)\dot{r}(x)-3\ddot{r}(x)^2}{4\dot{r}(x)^4}.\]

We discretize this equation with the pseudospectral method by expanding \(\phi_{n,l}(x)\) and \(\phi_{n,l}(x)/\dot{r}(x)\) as

\[\begin{split}\phi_{n,l}(x) &= \sum_{j=0}^N \phi_{n,l}(x_j) g_j(x), \\ \frac{\phi_{n,l}(x)}{\dot{r}(x)} &= \sum_{j=0}^N \frac{\phi_{n,l}(x_j)}{\dot{r}(x_j)} g_j(x).\end{split}\]

Inserting these expansions into Eq. symmetric_radial_TISE we have that

\[\sum_{j=0}^N \left(-\frac{1}{2} \frac{1}{\dot{r}(x)} \frac{\phi_{n,l}(x_j)}{\dot{r}(x_j)} g^{\prime \prime}_j(x) + V(r(x)) \phi_{n,l}(x_j) g_j(x) \right) = \epsilon_{n,l} \sum_{j=0}^N \phi_{n,l}(x_j) g_j(x)\]

Next, we multiply through with \(g_i(x)\) and integrate over \(x\),

\[\sum_{j=0}^N \left(-\frac{1}{2} \frac{\phi_{n,l}(x_j)}{\dot{r}(x_j)} \int \frac{g_i(x)}{\dot{r}(x)} g^{\prime \prime}_j(x) dx + \phi_{n,l}(x_j) \int g_i(x) V(r(x)) g_j(x) dx \right) = \epsilon_{n,l} \sum_{j=0}^N \phi_{n,l}(x_j) \int g_i(x) g_j(x) dx\]

The integrals are evaluated with by quadrature and using the property \(g_j(x_i) = \delta_{i,j}\) we have that

\[\begin{split}\int g_i(x) g_j(x) dx &= \sum_{m=0}^N g_i(x_m) g_j(x_m) w_m = \sum_{m=0} w_m \delta_{i, m} \delta_{j,m} = w_i \delta_{i,j}, \\ \int \frac{g_i(x)}{\dot{r}(x)} g^{\prime \prime}_j(x) dx &= \sum_{m=0} g_i(x_m) V(r(x_m)) g_j(x_m) w_m = w_i V(r(x_i)) \delta_{i,j}, \\ \int \frac{g_i(x)}{\dot{r}(x)} g^{\prime \prime}_j(x) dx & \underbrace{=}_{???} \sum_{m=1}^{N-1} \frac{g_i(x_m)}{\dot{r}(x_m)} g^{\prime \prime}_j(x_m) = w_i \frac{g^{\prime \prime}_j(x_i)}{\dot{r}(x_i)}, \ \ i=1,\cdots,N-1.\end{split}\]

Thus, for the interior grid points,

\[\sum_{j=1}^{N-1} \left(-\frac{1}{2} \frac{\phi_{n,l}(x_j)}{\dot{r}(x_j)} w_i \frac{g^{\prime \prime}_j(x_i)}{\dot{r}(x_i)} + \phi_{n,l}(x_j) w_i V(r(x_i)) \delta_{i,j} \right) = \epsilon_{n,l} \sum_{j=1}^{N-1} \phi_{n,l}(x_j) w_i \delta_{i,j}.\]

Using the expressions for the \(g_j^{\prime \prime}(x_i)\) we can write this as (notice that the weights \(w_i\) cancels)

\[\sum_{j=1}^{N-1} \left(-\frac{1}{2} \frac{\tilde{g}^{\prime \prime}_j(x_i) P_N(x_i)}{\dot{r}(x_i) \dot{r}(x_j)} \frac{\phi_{n,l}(x_j)} {P_N(x_j)} \right) + \phi_{n,l}(x_i) V(r(x_i)) = \epsilon_{n,l} \phi_{n,l}(x_i).\]

Furthermore, dividing through with \(P_N(x_i)\), we have that

\[\sum_{j=1}^{N-1} \left(-\frac{1}{2} \frac{\tilde{g}^{\prime \prime}_j(x_i)}{\dot{r}(x_i) \dot{r}(x_j)} \tilde{\phi}_{n,l}(x_j) \right) + V(r(x_i))\tilde{\phi}_{n,l}(x_i) = \epsilon_{n,l} \tilde{\phi}_{n,l}(x_i),\]

where we have defined

\[\tilde{\phi}_{n,l}(x_i) \equiv \frac{\phi_{n,l}(x_i)}{P_N(x_i)}.\]