Expectation value for $\frac{1}{r}$ in the $|nlm\rangle$ state of hydrogen I'm trying to derive $\left\langle \frac{1}{r} \right\rangle = \frac{1}{n^2 a_0}$ (where $a_0$ is the Bohr radius) for the $| nlm\rangle$ state of hydrogen.
I've separated the radial and angular parts of the hydrogen wavefunction and split up the integration to yield
\begin{align}\left\langle \frac{1}{r} \right\rangle &= \int R^*_{nl}(r) {Y^{m}_{l}}^*(\theta, \phi) \frac{1}{r} R_{nl}(r) Y_{l}^{m}(\theta, \phi) r^2 \sin \theta ~\mathrm dr ~\mathrm d\theta~\mathrm d\phi\\ &= \int {Y^{m}_{l}}^*(\theta, \phi) Y_{l}^{m}(\theta, \phi) \sin \theta ~\mathrm d\theta ~\mathrm d\phi \int \limits_{0}^{\infty} r R^*_{nl}(r) R_{nl}(r) ~\mathrm dr\\ &= \int \limits_{0}^{\infty} r R^*_{nl}(r) R_{nl}(r) ~\mathrm dr\end{align}
But from here, I'm not sure how to continue. Since $R_{nl}(r)$ is defined as
$$R_{nl}(r) = \left( \frac{2 Z}{n a_0} r \right)^l \sum \limits_{k = 0}^{n - l - 1} a_k \left( \frac{2 Z}{n a_0} r \right)^k e^{-Zr/n a_0}$$
where
$$a_{k + 1} = \frac{k + l + 1 - n}{(k + 1)(k + 2l + 2)}a_k$$
How do I deal with the sum inside the integral (especially a sum that requires recursion to compute coefficients)?
Have looked this source so far.
 A: That sum in your radial wavefunction is proportional to a generalized Laguerre polynomial $L_m^{(\alpha)}(x)$, and the more usual way to write the radial wavefunction is
$$R_{nl}(r)=N_{nl}\,\rho^l\,L_{n-l-1}^{(2l+1)}(\rho)\,e^{-\rho/2}$$
where
$$\rho\equiv\frac{2Z}{na_0}r$$
is a dimensionless radial coordinate and
$$N_{nl}\equiv\sqrt{\left(\frac{2Z}{na_0}\right)^3\frac{(n-l-1)!}{2n(n+l)!}}$$
is a wavefunction normalization factor that makes
$$\int_0^\infty r^2 R_{nl}^*(r)R_{nl}(r)dr=1.$$
The expectation value you want is thus
$$\begin{align}
\left\langle\frac1r\right\rangle&=\int_0^\infty rR_{nl}^*(r)R_{nl}(r)dr\\
&=N_{nl}^2\left(\frac{na_0}{2Z}\right)^2\int_0^\infty\rho^{2l+1}[L_{n-l-1}^{(2l+1)}(\rho)]^2e^{-\rho}d\rho\\
&=\frac{2Z}{na_0}\frac{(n-l-1)!}{2n(n+l)!}\int_0^\infty\rho^{2l+1}[L_{n-l-1}^{(2l+1)}(\rho)]^2e^{-\rho}d\rho\\
&=\frac{Z}{n^2a_0}\frac{(n-l-1)!}{(n+l)!}\int_0^\infty\rho^{2l+1}[L_{n-l-1}^{(2l+1)}(\rho)]^2e^{-\rho}d\rho.
\end{align}$$
Generalized Laguerre polynomials have been studied for about 150 years, and the integral above is very well known:
$$\int_0^\infty x^\alpha\,e^{-x}\,L_n^{(\alpha)}(x)\,L_m^{(\alpha)}(x)\,dx=\frac{(n+\alpha)!}{n!}\delta_{nm}.$$
Taking $n=m$ and substituting $\alpha\to 2l+1$ and $m\to n-l-1$, we find
$$\begin{align}
\left\langle\frac1r\right\rangle&=\frac{Z}{n^2a_0}\frac{(n-l-1)!}{(n+l)!}\int_0^\infty\rho^{2l+1}[L_{n-l-1}^{(2l+1)}(\rho)]^2e^{-\rho}d\rho\\
&=\frac{Z}{n^2a_0}\frac{(n-l-1)!}{(n+l)!}\frac{(n+l)!}{(n-l-1)!}\\
&=\frac{Z}{n^2a_0}.
\end{align}$$
A: While you should learn to compute the integrals, let me point out in this additional answer that one can compute this without evaluating any integrals. 
The Hamiltonian is:
$$H = \frac{p^2}{2m} + \frac{e^2}{r}$$
the energy eigenvalue of the $\left|n,l,m\right>$ eigenstate is
$$E_{nlm} = \frac{m e^4}{2\hbar^2 n^2}\tag{1}$$
here I've written this explicitly in terms of the constants that appear in the Hamiltonian. If we add the term $$V = \frac{g}{r}$$ as a perturbation to the Hamiltonian, then we know from first order perturbation theory that the perturbed energy eigenvalues are given to first order in $g$ by:
$$E(g) = E_{nlm} + g \left<nlm\right|\frac{1}{r}\left|nlm\right>+\mathcal{O}(g^2)\tag{2}$$
We can also compute the perturbed energy directly by modifying the charge $e$. We can write:
$$H + V =  \frac{p^2}{2m} + \frac{e'^2}{r}$$
with:
$$e'^2 = e^2 + g$$
From Eq. (1) it follows that:
$$E(g) = \frac{me'^4}{2\hbar^2} = E_{nlm} + g\frac{me^2}{\hbar^2 n^2} + \mathcal{O}(g^2)$$
Comparing this to Eq. (2) yields:
$$\left<nlm\right|\frac{1}{r}\left|nlm\right> = \frac{me^2}{\hbar^2 n^2}=\frac{1}{a_0 n^2}$$
A: Not a full answer but I think you can switch  sigma and integral signs  around. 
Also look at Sums and Integrals
That is:
$${\displaystyle \int \sum _{r=a}^{b}f\left(r,x\right)\,dx=\sum _{r=a}^{b}\int f\left(r,x\right)\,dx}$$
This is simply because:
\begin{align}
 \int \sum _{r=a}^{b}f(r,x)\,dx &=\int f\left(a,x\right)+f((a+1),x)+f((a+2),x)+\dots  +f((b-1),x)+f(b,x)\,dx\\
 &=\int f(a,x)\,dx+\int f((a+1),x)\,dx+\dots +\int f((b-1),x)\,dx+\int f(b,x)\,dx\\ &=\sum _{r=a}^{b}\int f(r,x)\,dx
\end{align}
