Solution (source code)

= Solution

A pure <stationary state> of a time-independent <Hamiltonian operator> is $\Psi(\mathbf r,t)=\psi(\mathbf r)e^{-iEt/\hbar}$. Its <probability density> and all expectations of time-independent <observables> are constant. This does not require zero <probability current>.

For real scalar <potential energy> and no <vector potential>, substitute $i\hbar\dot\Psi=-(\hbar^2/2m)\nabla^2\Psi+V\Psi$ and its <complex conjugate> into $\dot\rho=\Psi^*\dot\Psi+\dot\Psi^*\Psi$. The potential terms cancel:
$$
\dot\rho=\frac{i\hbar}{2m}\big(\Psi^*\nabla^2\Psi-\Psi\nabla^2\Psi^*\big)
=-\nabla\cdot\left[\frac{\hbar}{2im}(\Psi^*\nabla\Psi-\Psi\nabla\Psi^*)\right].
$$
Thus the <probability current> is $\mathbf j=(\hbar/m)\operatorname{Im}(\Psi^*\nabla\Psi)$. In a <stationary state> $\nabla\cdot\mathbf j=0$: flow may circulate without changing the <density>.

Here only the azimuthal <phase> contributes an imaginary logarithmic <derivative>. Since $\nabla\phi=\hat{\boldsymbol\phi}/(r\sin\theta)$,
$$
\boxed{\mathbf j=\frac{\hbar|\Psi|^2}{mr\sin\theta}\hat{\boldsymbol\phi}
=\frac{\hbar A^2}{m}r e^{-r/a_0}\sin\theta\,\hat{\boldsymbol\phi}}.
$$
There is no radial or polar flow; the azimuthal current is independent of $\phi$, so its <divergence> vanishes. The apparent coordinate singularity at $\sin\theta=0$ cancels in the final expression.

The z-component of <orbital angular momentum> follows from the moment of the <mass> current:
$$
\langle L_z\rangle
=\int mr\sin\theta\,j_\phi\,dV
=\hbar A^2
\left(\int_0^\infty r^4e^{-r/a_0}dr\right)
\left(\int_0^\pi\sin^3\theta\,d\theta\right)
\left(\int_0^{2\pi}d\phi\right)
=\boxed{\hbar}.
$$
Indeed $A^2(24a_0^5)(4/3)(2\pi)=1$. This agrees with $L_z=-i\hbar\partial_\phi$ acting on the <wavefunction>.