I'm interested in the general formulas that give the exact uncertainties $\Delta r$ and $\Delta p_r$ (the radial momentum) for all stationary states $|n,l, m \rangle$ (or $\psi_{nlm}(r, \theta, \varphi)$), as functions of the quantum numbers $n = 1, 2, 3, \ldots, \infty$ and $l = 0, 1, 2, \ldots, n - 1$. The radial momentum hermitian operator I'm considering here is $$ p_r = -\, i \hbar \, \frac{1}{r} \, \frac{\partial }{\partial r} (\, r \quad ) = -\, i \hbar \Bigl( \frac{\partial }{\partial r} + \frac{1}{r} \Bigr). \tag{1} $$ This operator is hermitian: $p_r^{\dagger} = p_r$, and satisfies the usual commutator: $[r, p_r] = i \hbar$, but may have a physical interpretation issue (see the comments below).
I'm able to calculate the uncertainties in the special case of $l = n - 1$ (since the Laguerre polynomials are just $L_0^k(u) = 1$ in this case): \begin{align} \Delta r &= \frac{n a}{2 \mathrm{Z}} \, \sqrt{2 n + 1}, \tag{2} \\[2ex] \Delta p_r &= \frac{\mathrm{Z} \hbar}{n a} \, \frac{1}{\sqrt{2 n - 1}}. \tag{3} \end{align} The Heisenberg relation is thus satisfied (of course): $$\tag{4} \Delta r \, \Delta p_r = \frac{\hbar}{2} \, \sqrt{\frac{2 n + 1}{2 n - 1}} \ge \frac{\hbar}{2}. $$ The equality is realized for $n \gg 1$. But then, these expressions are only valid for $l = n - 1$. What about the general case?
Are the stationary states position and momentum uncertainties known as functions of $n$ and $l$?
It's easy to calculate $\langle \, p_r \rangle$ for all stationary states $\psi_{n l m}(r, \theta, \varphi)$: $$\tag{5} \langle \, p_r \rangle = 0. $$ But for $\langle \, p_r^2 \rangle$, it's much harder because of the Laguerre polynomials. I can only show this integral (from an integration by parts. $\phi_{n l}(r)$ is the radial part of the atomic wave function): $$\tag{6} \langle \, p_r^2 \rangle = \hbar^2 \int_0^{\infty} \Bigl( \frac{d\, }{d r} \bigl( r \, \phi_{n l}(r) \bigr) \Bigr)^{2} \, dr. $$
The comments below are objecting that the radial momentum (1) isn't a good observable because of the spectral decomposition issue. But then, the operator $p_r^2$ is showing in the atomic Hamiltonian, which is a good Hermitian operator associated to a well behaved observable. So while $p_r$ may be problematic as an observable, $p_r^2$ isn't. What I'm interested in is the average (6), so this is a perfectly valid question. I guess the simplest way to get (6) is to find $\langle \, r^{-1} \rangle$ and $\langle \, r^{-2} \rangle$ (for a general stationary state $|n, l, m \rangle$) and to substract their contributions to the energy levels: $$\tag{7} \frac{1}{2 m} \langle \, p_r^2 \rangle = E_n - \frac{\hbar^2 \, l (l + 1)}{2 m} \langle \, r^{-2} \rangle + k \mathrm{Z} e^2 \langle \, r^{-1} \rangle. $$ So the question remains the same and I think the objections below are irrelevant.