In today's exercises, we complete the tasks from yesterday, which unfortunately had ambiguous and confusing instructions!
Here, our goal is to do some simple exercises to help us build useful tools such as subroutines to calculate HO wf's and matrix elements $\langle n'l|V|nl\rangle$ that are useful building blocks for more complicated calculations (e.g., Hartree-Fock calculations of many-body systems). Our ultimate test will be to see if we can get the ground state energy of the hydrogen atom (-.5 in atomic units) when we solve the Schr\“odinger equation by matrix diagonalization in the HO basis.
In the code
here, find the subroutine that computes the generalized laguerre polynomials $L_n^{l+1/2}(x)$, “laguerre_general”. NOTE: this is not a standalone code. Rather, the purpose is to see the simple recursion algorithm for calculating the Laguerre polynomials, and by extension $R_{nl}$, which you can then implement according to your own taste.
Recall that the HO wavefunctions are given by
\begin{equation}
R_{nl}(r) = \frac{A_{nl}}{b^{3/2}}\xi^le^{-\xi^2/2}L_{n}^{l+1/2}(\xi^2)\,,
\end{equation}
where the oscillator length $b\equiv \sqrt{\hbar/(m\omega)}$, and the dimensionless variable $\xi\equiv r/b$. The normalization constant is given by
\begin{equation}
A_{nl}=\sqrt{\frac{2^{n+l+2}n!}{\pi^{1/2}(2n+2l+1)!!}}
\end{equation}
Using the uploaded code for guidance, make your own subroutine to calculate $L_{n}^{\alpha}(x)$. Then, using the definition of the HO wf's, make a subroutine/function to calculate $R_{nl}(r)$.
\begin{equation}
\int_{0}^{\infty}R_{nl}(r)R_{n'l}(r)r^2dr= \delta_{nn'}
\end{equation}
As discussed in the HO notes, the “best” type of Gaussian quadrature over HO wf's is Gauss-Laguerre. However, as you see in the code
here, you can use Gauss-Legendre as well. You need to play around to make sure your upper limit $r_{max}$ of your discretization range is sufficiently large. (Note: there are plenty of freely downloadable Gauss-Legendre subroutines. Please feel free to use one instead of writing it from scratch.) In this way, the orthogonality check takes the form
\begin{equation}
\int_{0}^{\infty}R_{nl}(r)R_{n'l}(r)r^2dr\approx \sum_{i=1}^{N}w_i r_i^2 R_{nl}(r_i)R_{n'l}(r_i)\,,
\end{equation}
where the sum is over the N mesh points and $w_{i}$ are the quadrature weights.
To test your $\langle nl|V|n'l\rangle$ routine, apply it to the coulomb potential in the hydrogen atom problem. Use atomic units ($\hbar=e=1/4\pi\epsilon_0=m_e=1$). Now, build the Hamiltonian matrix $\langle nl|H|n'l\rangle$ (you may take $l=0$) for the hydrogen atom, keeping states $n,n'< N_{max}$. Use the analytic expressions for the K.E. matrix elements in
ho_spherical.pdf. Now diagonalize the matrix to get the ground state energy.
For a given $N_{max}$, repeat the calculation for many different values of oscillator length and plot the ground state versus $b$. Repeat for larger values of $N_{max}$ and put the $E$ versus $b$ curves on the same plot. Are they behaving according to expectation? (Hint: the diagonalization in a finite basis is variational.) Should results depend on $b$ as $N_{max}\rightarrow\infty$?