Quantization of a particle on a spherical surface Suppose we have a particle of mass $m$ confined to the surface of a sphere of radius $R$. The classical Lagrangian of the system is 
$$L = \frac{1}{2}mR^2 \dot{\theta}^2 + \frac{1}{2}m R^2 \sin^2 \theta \dot{\phi}^2 $$ 
The canonical momenta are 
$$P_\theta = \frac{\partial L }{\partial \dot{\theta }} = m R^2 \dot{\theta }$$ and $$P_\phi = \frac{\partial L }{\partial \dot{\phi }} = m R^2 \sin^2 \theta \dot{\phi }$$
The Hamiltonian is 
$$H = \frac{P_\theta^2}{2 m R^2} + \frac{P_\phi^2}{2 m R^2 \sin^2\theta }$$
Now start to quantize the system. We replace $P_\theta $ and $P_\phi $ as $-i\hbar \frac{\partial}{\partial \theta}$ and $-i\hbar \frac{\partial}{\partial \phi} $, respectiely, to obtain 
$$H = -\frac{ \hbar^2 \partial^2}{2 m R^2 \partial \theta^2} - \frac{\hbar^2 \partial^2 }{2 m R^2 \sin^2\theta \partial \phi^2 } $$ 
This is apparently wrong, it should be the total angular momentum! 
So what is the right procedure to quantize a system, especially a system in curvilinear coordinates?
 A: Conventional wisdom (as stated in the textbooks of Shankar, or Griffiths, for example) says to avoid quantizing operators in curvilinear coordinates whenever possible. Better to quantize the cartesian operators $p_x, p_y, p_z$ and then switch to curvilinear coordinates in the quantum theory. Refer to this discussion by Professor Robert Jaffe at MIT:[1]

The canonical quantization method becomes complicated and subtle when
  one tries to apply it to coordinate systems that include singular
  points. A familiar example is spherical polar coordinates
  $(r,\theta,\phi)$. The origin, $r=0$, is a singular point for
  spherical polar coordinates---for example, $\theta$ and $\phi$ are not
  defined at $r=0$. If you follow the canonical formalism through from
  Lagrangian to canonical momenta $(p_r,p_{\theta},p_{\phi})$ to
  Hamiltonian, to canonical commutators, a host of difficulties arise.
  Although it is possible to sort them out by insisting that all the
  canonical momenta be Hermitian operators, it is considerably easier to
  quantize the system in Cartesian coordinates and make the change to
  spherical polar coordinates at the quantum level. This is the path
  taken in most elementary treatments of quantum mechanics in three
  dimensions: the operator $p^2 = p_1^2 + p_2^2 +p_3^2$ is recognized as
  the Laplacian in the coordinate representation $(p_j \to -i\hbar
> \partial/\partial x_j \implies p^2 \to \nabla^2)$ and the
  transformation to polar coordinates is made by writing the Laplacian
  and the wavefunction in terms of $r$, $\theta$, and $\phi$. As a rule
  of thumb, the canonical approach becomes cumbersome when the classical
  coordinates and/or momenta do not range over the full interval from
  $-\infty$ to $+\infty$.

[1] R.L. Jaffe, "Canonical Quantization and Application to
the Quantum Mechanics of a Charged
Particle in a Magnetic Field", Supplementary Notes for MIT’s Quantum
Theory Sequence, February 2007
A: In a nutshell, the problem with OP's choice of operators $\hat{p}_j$ and $\hat{H}$ is that they are not selfadjoint wrt. to the pertinent measure $\mu$. In other words, the usual integration by parts method to prove selfadjointness does not work. 
Here are some more details. Let us put the constants $m=1=R$ for simplicity. Then OP's Lagrangian becomes of the form
$$\tag{1} L~=~\frac{1}{2}g_{ij}~\dot{x}^i\dot{x}^j,$$
with coordinates $x^1\equiv\theta$, $x^2\equiv\phi$, and metric tensor
$$\tag{2}  g_{ij}~=~ \begin{pmatrix} 1 & 0 \\ 0 &\sin^2\theta \end{pmatrix}. $$
Classically, the Lagrangian momenta are 
$$\tag{3} p_i ~=~g_{ij}~\dot{x}^j,$$ 
and the Hamiltonian is
$$\tag{4} H~=~\frac{1}{2} g^{ij} ~p_i p_j. $$
The volume form in configuration space is 
$$\tag{5} \mu ~=~ \sqrt{g}~ \mathrm{d}x^1 \wedge \mathrm{d}x^2  
~=~ \sin{\theta} ~\mathrm{d}\theta \wedge \mathrm{d}\phi. $$ 
The Hilbert space is $L^2(S^2,\mu)$. What is the Schrödinger representation of the momentum operators? Well, now we run into operator ordering ambiguities. The momentum operators should as a minimum satisfy (i) the CCR, and (ii) be selfadjoint wrt. to the measure (5). One idea to ensure this is to use
$$\tag{6} \hat{p}_j~=~ \frac{\hbar}{i\sqrt[4]{g}} \frac{\partial}{\partial x^j} \sqrt[4]{g}. $$
Similarly, we can choose a selfadjoint Hamiltonian operator to be the Laplace-Beltrami operator:
$$\tag{7}  \hat{H}~=~-\frac{\hbar^2}{2}\Delta
~=~ -\frac{\hbar^2}{2\sqrt{g}}\frac{\partial}{\partial x^i}\sqrt{g}~ g^{ij} \frac{\partial}{\partial x^j}~=~ \frac{1}{2\sqrt[4]{g}} \hat{p}_i\sqrt{g}~ g^{ij} ~\hat{p}_j\frac{1}{\sqrt[4]{g}} .
$$
In case of the two-sphere $S^2$, this Hamiltonian operator (7) leads to the square of the angular momentum $\hat{\bf L}^2$. Classically, the operators (6) and (7) reduce to the functions (3) and (4), respectively.
References:


*

*Bryce DeWitt, Supermanifolds, Cambridge Univ. Press, 1992; Section 6.7.

