Since the Laguerre polynomial is deeply connected to the hydrogen-like electron orbital. The radial solution is

\displaystyle R_{nl}(r) = A \frac{1}{r} \exp(- \frac{x}{2}) x^{l+1}L_{n-l-1}^{2l+1}(x)

\displaystyle x = \frac{2Z}{na_0} r

From the normalization condition, we can get the coefficient A

\displaystyle \int_0^{\infty} R_{nl}^2(r) r^2 dr

Beside of calculating the normalization factor, the general expectation value also need to evaluation of the integration of the Laguerre polynomial.

In here, we will show the calculation for 3 integrals:

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx,  n,\alpha\geq0

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx

Using the Rodrigues formula

\displaystyle L_n^\alpha(x) = \frac{1}{n!} x^{-\alpha} e^x \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = \frac{1}{n!} \int_0^{\infty} x^{m} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx

For m=0,

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx \\ = \frac{1}{n!} \int_0^{\infty} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx \\ = \frac{1}{n!} (F_{n-1}(\infty)-F_{n-1}(0) )

\displaystyle F_{n-1}(x) = \int_0^{\infty} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx \\ = \sum_{i=0}^{n-1} (-1)^i \begin{pmatrix} n-1 \\ i \end{pmatrix} \frac{(n+\alpha)!}{(\alpha+1+i)!} x^{\alpha+1+i} e^{-x}

It is obvious that F(\infty) = 0 due to the e^{-x}. Since \alpha>0, then x^{\alpha+1+i} = 0 at x = 0. Thus,

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = 0 ,  m=0

Using integration by path,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = [ x^{m} \frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) ]_0^{\infty} - (m) \int_0^{\infty} x^{m-1} \frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha})dx

For the same reason, the first term is zero for \alpha + m \geq 0. Repeat the integration by path k times,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^k \frac{(m)!}{(m-k)!} \int_0^{\infty} x^{m-k} \frac{d^{n-k}}{dx^{n-k}}(e^{-x}x^{n+\alpha})dx

When m \geq n || 0 > m, k can go to n, then,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^{n} \frac{(m)!}{(m-n)!} \int_0^{\infty} x^{m-n}(e^{-x}x^{n+\alpha})dx \\ = (-1)^{n} \frac{(m)!}{(m-n)!} (\alpha+m)!

When n > m > 0, k can only go to m. then,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^{m} (m)! \int_0^{\infty} \frac{d^{n-k}}{dx^{n-k}}(e^{-x}x^{n+\alpha})dx \\ = (-1)^{m} (m)!  \frac{d^{n-m-1}}{dx^{n-m-1}}(e^{-x}x^{n+\alpha}) = 0

For the similar result as F_{n-1}(x), the result is zero.

To summarize, the following formula suitable for all m > n || 0>m

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = (-1)^n (\alpha+m)! \begin{pmatrix} m \\ n \end{pmatrix}

To evaluate the integral,

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx

We need the know that the Laguerre polynomial can be expressed as,

\displaystyle L_n^\alpha(x) = \sum_{i=0}^{n} (-1)^i \begin{pmatrix} n+\alpha \\ n-i \end{pmatrix} \frac{x^i}{i!}

Thus, we can evaluate

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} x^i L_n^\alpha(x) dx

Then combine everything,

\displaystyle \int_0^{\infty} x^{\alpha+k+i} e^{-x} L_n^\alpha(x) dx = (-1)^n (\alpha+i+k)! \begin{pmatrix} i+k \\ n \end{pmatrix}

put inside the series expression,

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx = \sum_{i=0}^{n} (-1)^{i+n} \begin{pmatrix} n+\alpha \\ n-i \end{pmatrix} \begin{pmatrix} i+k \\ n \end{pmatrix} \frac{(\alpha+k+i)!}{i!}

In general, i + k \geq n || 0 > i+k .

For k = 0 , the sum only has 1 term for i = n

\displaystyle \int_0^{\infty} x^{\alpha} e^{-x} (L_n^\alpha(x))^2 dx = \frac{(\alpha+n)!}{n!}

Using the formula, for k=1, we have

\displaystyle \int_0^{\infty} x^{\alpha+1} e^{-x} (L_n^\alpha(x))^2 dx = \frac{(\alpha+n)!}{n!} (2n+\alpha +1)

To evaluate the expectation value of radial function

\displaystyle \langle R_{nl} | \frac{1}{r^m} | R_{nl} \rangle

We have to calculate the integral,

\displaystyle \int_0^\infty e^{-x} x^{2l+2-m} \left( L_{n-l-1}^{2l+1} (x) \right)^2 dx

For m = 0, we get the normalization constant,

\displaystyle A^2 \frac{na_0}{2Z} \int_0^{\infty} e^{-x} x^{2l+2} \left( L_{n-l-1}^{2l+1}(x) \right)^2 dx = 1

\displaystyle \langle \frac{1}{r} \rangle = \frac{1}{n^2} \frac{Z}{a_0}

\displaystyle \langle \frac{1}{r^2} \rangle = \frac{2}{n^3(2l+1)} \frac{Z^2}{a_0^2}

In general,

\displaystyle \langle \frac{1}{r^m} \rangle = \left(\frac{Z}{a_0}\right)^m \frac{2^{m-1}(n-l-1)!}{n^{m+1}(n+l)!} G_{nlm}

\displaystyle G_{nlm} = \sum_{i=0}^{n-l-1} (-1)^{n-l-1+i} \begin{pmatrix} n+l \\ n-l-1+i \end{pmatrix} \begin{pmatrix} 1-m+i \\ n-l-1 \end{pmatrix} \frac{(2l+2-m+i)!}{i!}

For the integral, for m \neq n

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx

When k = 0, when m < n , the integral is zero.

When k \neq 0, m<n

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx \\ = \sum_{i=0}^{m} \frac{(-1)^i}{i!}\begin{pmatrix} m+\beta \\ m-i \end{pmatrix} \int_0^{\infty} x^{\alpha+k+i} e^{-x} L_n^\alpha(x) dx \\ = \sum_{i=0}^{m} \frac{(-1)^{i+n}}{i!}\begin{pmatrix} m+\beta \\ m-i \end{pmatrix} \begin{pmatrix} k+i \\ n \end{pmatrix} (\alpha+k+i)!

The sum only non-zero when m\geq i > n-k or -k > i \geq 0

This shows that the Laguerre polynomial is orthogonal with weighting e^{-x} x^{\alpha}.