Hamiltonian operator in polar coordinates with momentum operators The Hamiltonian operator for a free non-relativistic particle looks like
$$
\hat{H}  = \frac{\hat{p}^2}{2m}   = -\frac{\hbar^2}{2m} \nabla^2. $$
In polar coordinates, the Laplacian expands to
$$
\hat{H}  = -\frac{\hbar^2}{2m} \left (\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} \right).  $$
The radial and angular momentum operators are 
$$
\begin{align*}
\hat{p}_r = \frac{\hbar}{i} \left(\frac{\partial}{\partial r} + \frac{1}{2r}\right) && \hat{p}_\theta = \frac{\hbar}{i}\frac{1}{r}\frac{\partial}{\partial \theta}.
\end{align*}$$
After squaring, summing and comparing to the Hamiltonian, we find that
$$ \hat{H} = -\frac{\hbar^2}{2m} ( \hat{p}_r^2 + \hat{p}_\theta^2 ) - \frac{\hbar^2}{8mr^2}. $$
In classical mechanics, we expect that $ p^2 = p_r^2 + p_\theta^2 $, does this not hold in quantum mechanics? Why does this curious $- \frac{\hbar^2}{8mr^2}$ potential appear, does it have any significance?
Addendum
To clarify the choice of radial momentum operator, consider the naive $ \hat{p}'_r = \frac{\hbar}{i} \frac{\partial}{\partial r} $. Taking the adjoint we find that
$$ \langle \phi, \hat{p}'_r \psi \rangle = \int d\theta \int_0^{+\infty} r dr \phi^* \frac{\hbar}{i} \frac{\partial}{\partial r} \psi = - \int d\theta \int_0^{+\infty} r dr \frac{\hbar}{i} \left( \frac{\partial \phi^*}{\partial r} + \frac{\phi^*}{r} \right)\psi = \int d\theta \int_0^{+\infty} rdr \frac{\hbar}{-i} \left( \frac{\partial}{\partial r} + \frac{1}{r} \right) \phi^* \psi \neq \langle \hat{p}'_r \phi, \psi \rangle$$
With an extra term $\frac{\hbar}{i} \frac{1}{2r}$ it is self-adjoint. This term is different from the one in spherical coordinates by a factor of 2.
 A: Coordinate transforms do not work the same in quantum mechanics as in classical mechanics. Specifically, canonical quantization is not invariant with respect to most transformations in phase space, or even with respect to just spatial coordinate transformations. At the end of the day, a simple recipe is to just transform the Schrödinger equation in Cartesian coordinates as a partial differential equation without any interpretation of the new terms as momenta with respect to the new coordinates. Let me illustrate my point below.

When discussing coordinate transformations in QM, you can take one of the two approaches:


*

*You do quantization in Cartesian coordinates $x^i$, and are transforming the Schrödinger equation for $\psi^{(x)}(x^i)$ to a new set of coordinates $q^i$. That is, you carefully follow the transformation of the differential operator $\sum_i \partial^2/\partial {x^i}^2$ and the volume element $d^3 x$ according to the usual mathematical rules. As mentioned, there is not really an interpretation in terms of canonical momenta of the new equations, but the nice thing is $|\psi^{(x)}(q^i)|^2$ has the meaning of a probability density per actual physical volume in any coordinates. 

*You make a coordinate transform on the coordinate part of the phase space to some new coordinates $q^i$. Then you canonically quantize, which replaces any of your canonical momenta only as $p_i \to \hat{P}_i = -i \hbar \partial/\partial q^i$ in the Hamiltonian. The resulting Schrödinger equation holds for a wave function $\psi^{(q)}(q^i)$ that corresponds to a coordinate volume $d^3 q$.


Consider the example of the spherical polar coordinates. Approach 1. will give the Schrödinger equation as
$$-i\hbar \frac{\partial \psi^{(x)}}{\partial t} = -\frac{\hbar^2}{2m} \left (\frac{1}{r^2} \frac{\partial}{\partial r^2} \left( r^2 \frac{\partial}{\partial r} \right) + \frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta} \left(\sin \theta \frac{\partial}{\partial \theta} \right) + \frac{1}{r^2 \sin^2 \! \theta} \frac{\partial^2}{\partial \varphi^2}\right)\psi^{(x)} + V(r,\theta, \varphi)\psi^{(x)}$$
The quantity $|\psi^{(x)}(r,\theta,\varphi)|^2$ is here the density per physical volume $r^2 \sin\theta\, dr d\theta d \varphi$.
For approach 2. we first need to transform the classical Hamiltonian to polar coordinates
$$ H = \frac{1}{2m}\left(p_r^2 + \frac{p_\theta^2}{r^2} + \frac{p_\varphi^2}{r^2 \sin^2\! \theta}\right) + V(r,\theta,\varphi) $$
Now we have $\hat{P}_r = -i\hbar \partial/\partial r, \hat{P}_\theta = -i\hbar \partial/\partial \theta, \hat{P}_\varphi = -i\hbar \partial/\partial \varphi$, and the Schrödinger equation obviously is
$$-i\hbar \frac{\partial \psi^{\rm (sph)}}{\partial t} = \hat{H}\psi = -\frac{\hbar^2}{2m} \left (\frac{\partial^2}{\partial r^2}  + \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} + \frac{1}{r^2 \sin^2 \! \theta} \frac{\partial^2}{\partial \varphi^2}\right)\psi^{\rm (sph)} + V(r,\theta, \varphi)\psi^{\rm (sph)}$$
Note that now the meaning of $|\psi^{\rm (sph)}|^2$ is that of a density per coordinate volume $dr d\varphi d\vartheta$. This would suggest that the two wavefunctions are related by a factor $r \sqrt{\sin \theta}$ and perhaps a phase factor. However, if you try to transform one equation into another in this way, you see that they are simply inequivalent. 
Only one of the two equations can be true, because they give different experimental predictions. It turns out that the correct one is the equation obtained by the Cartesian approach 1. More on the question of coordinates for quantization can be found in this answer.
A: The main reason you can't square the components here directly is that the basis vectors $\hat{r}$ and $\hat{\theta}$ are changing with the coordinates. The second reason is that the metric of this coordinate is not $\mathbf{1}$. (This coordinate is a curved orthogonal coordinate but cartesian coordinate is a flat coordinate.) Thus, when taking the derivative of a vector field $\mathbf{F}$, you have $$\partial_j\mathbf{F}=\partial_j F^i\mathbf{g}_i+F^i\partial_j\mathbf{g}_i=\left(\partial_jF^i+F^k\Gamma^i_{jk}\right)\mathbf{g}_i$$ Where the additional second term originates from the difference of basis vectors at different points.
Thus $\nabla\cdot\mathbf{F}=\partial_iF^i=\left(\partial_iF^i+F^k\Gamma^i_{ik}\right)\mathbf{g}_i=\frac{1}{\sqrt{g}}\partial_k(\sqrt{g}F^k)$. The $\nabla^2$ operator is then $\nabla^2=\nabla\cdot\nabla=\partial_i\partial^i=\frac{1}{\sqrt{g}}\partial_k(\sqrt{g}\partial^k)=\frac{1}{\Pi_kA_k}\partial_i(\frac{\Pi_jA_j}{A_i^2}\partial_i)$. Where $A_i$ are the Lame constants.
Here the $A_i$ are $1$ and $r$ respectively, simply plug into the formula for $\nabla\cdot$ and $\nabla^2$ You get $\nabla_r=\frac{1}{r}\partial_r r=\frac{1}{r}+\partial_r$ (you had a coefficient wrong) and $\nabla_\theta=\frac{1}{r}\partial_\theta$. And $$\nabla^2=\frac{1}{r}\left(\partial_r \left(r\partial_r\right)+\partial^2_\theta\right)$$
