I want to compute quickly (using maybe some scaling arguments) $\langle n^{'}l^{'}m^{'}|r^k|nlm\rangle$, where $k \in I$. $|nlm \rangle$ is the eigenfunction of the Hydrogen atom ($H$).
Example:
Quickly compute $\langle r\rangle_{nlm}$ and the likes.
Some quick help would be highly appreciated.
This looks to be a messy calculation for the most general case, and a direct attempt based on special-functions properties of the Laguerre wavefunctions is likely to simply falter and die in the not-quite-right forms of the integrals.
These matrix elements are calculated in terms of recursion relations in
Matrix-element calculations for hydrogenlike atoms. M.L. Sánchez, B. Moreno, and A. López Piñeiro. Phys. Rev. A 46 no. 11, 6908 (1992)
and the authors provide software to calculate them in
HYDMATEL: a code to calculate matrix elements for hydrogen-like atoms. M.L. Sánchez and A.López Piñeiro. Comput. Phys. Commun. 75 no. 1-2, 185 (1993)
with the software available at
HYDMATEL: a code to calculate matrix elements for hydrogen-like atoms. Computer Physics Communications Program Library, Id ACLO v1.0
under the CPC standard license.
This answer gives an analytical approach for diagonal matrix elements.
First of all, since $r^k$ is spherically symmetric you can immediately integrate the angular parts: \begin{equation} \left< n' l' m' | r^k | n l m \right> = \left< n'l \right\| r^k \left\| nl \right> \delta_{l',l} \delta_{m',m}, \end{equation} where $\left< n'l \right\| r^k \left\| nl \right>$ is called the reduced matrix element.
Diagonal matrix elements: expectation values of $r^k$
Diagonal matrix elements $\langle r^k \rangle = \left< nl \right\| r^k \left\| nl \right>$ can be calculated with a combination of the Hellmann-Feynman theorem and Kramers' relation without the need to perform any integral explicitly.
The Hellmann-Feynman theorem is given by \begin{equation} \frac{dE}{d\lambda} = \left< \psi(\lambda) \right| \frac{\partial H}{\partial \lambda} \left| \psi(\lambda) \right>, \end{equation} where $\left|\psi(\lambda) \right>$ is a normalized ket for all $\lambda$. In this case you take $\lambda = l$ (note that $n=n(l)$) for $k=-2$ and $\lambda = Z$ for $k=-1$. Here $l$ is the angular momentum quantum number and $Z$ the atomic number.
Kramers' relation is a recursion relation from which all other $\langle r^k \rangle$ for integer $k$ can be found. It can be derived (see link) from the radial Schroedinger equation of the hydrogen-like atom and is given by \begin{equation} \frac{k}{4} \left[ \left( 2l + 1 \right)^2 - k^2 \right] a_0^2 \langle r^{k-2} \rangle - Z \left( 2k + 1 \right) a_0 \langle r^{k-1} \rangle + Z^2 \frac{k+1}{n^2} \langle r^k \rangle = 0, \end{equation} where $a_0$ is the Bohr radius. Note that this relation by itself gives all $\langle r^k \rangle$ for $k>-2$ since $\langle r^0 \rangle = 1$.
Off-diagonal matrix elements
For the off-diagonal matrix elements, I looked in the 8th edition of Gradshteyn and Ryzhik (page 817) for an analytical expression of the relevant integrals, but none is given there. Likely, there are no analytical expressions for the off-diagonal matrix elements. This is also confirmed by the paper linked in the answer of Emilio Pisanty.
Given that the eigenstates of the hydrogen atom can be separated into a tensor product of angular momentum eigenstates and a radial solution, the Wigner-Eckart theorem helps to calculate the matrix element of any tensor-like operator (a scalar as $r^k$ is a trivial case) between states of the form $|lm\rangle$. The rest is an integral over the radial components that can be calculated using the fact that the Laguerre polynomials (contributing to the radial component) form an orthonormal basis in $L^2$. Obviously, some integrals may be complicated to work out, but that is a different matter.
