Consider a particle constrained to move within a fixed Riemann manifold $(M, g)$ and minimally coupled to a fixed 1-form field $A$ and a fixed scalar field $U$. The canonical Hamilton function $H_C$ for this system is given in local coordinates by
$$H_C(q, p, t) = \frac{1}{2m}g^{ij}(q)\big(p_i - A_i(q, t)\big)\big(p_j - A_j(q, t)\big) + U(q, t).$$
My question concerns the quantization of this system. Conventional wisdom dictates that we forgo the use of non-Cartesian coordinates for as long as possible when quantizing a mechanical system. However, that does not seem possible here in general. Let us try a different approach, namely the phase space formulation. We replace operator multiplication with the star product of phase space functions:
$$f \star g = \exp\bigg(i\frac{\hbar}{2}P\bigg)(f,g)$$
where the bidifferential operator $P$ is given in local coordinates by
$$P(f,g) = \frac{\partial f}{\partial q^i}\frac{\partial g}{\partial p_i} - \frac{\partial f}{\partial p^i}\frac{\partial g}{\partial q_i}.$$
We replace operator commutation with the Moyal bracket of phase space functions:
$$\big\{\{f, g\}\big\} = -i\frac{1}{\hbar}(f \star g - g \star f).$$
The state of the system is defined by the Wigner function $W$, a quasiprobability density on phase space, whose time evolution is given by a quantum modification of Liouville flow:
$$\frac{\partial W}{\partial t} = -\big\{\{W, H\}\big\}.$$
My uncertainty centers around the Hamilton function $H$ that appears in the time evolution above. According to the Wikipedia article on the phase space formulation (https://en.wikipedia.org/wiki/Phase-space_formulation), this $H$ is "most often identical to the original Hamiltonian", which suggests that we choose $H = H_C$. I have my doubts about this, though, considering that there may be operator-ordering ambiguities on the Hilbert space side whose nature might suggest an altogether different choice for $H$.
What, then, ought we choose for $H$? I understand there may not be a clear answer, since the classical/quantum correspondence is not unique. In that case, have any experiments shed light on what is the most natural choice?
Edit: Specialize to the case of Euclidean space in Cartesian coordinates. The canonical Hamilton function becomes
$$H_C(x, p, t) = \frac{1}{2m}\delta^{ij}\big(p_i - A_i(x, t)\big)\big(p_j - A_j(x, t)\big) + U(x, t)$$
Expanding this out, we have
$$H_C = \frac{1}{2m}\delta^{ij}\big(p_ip_j - A_i(x, t)p_j - p_iA_j(x, t) + A_i(x, t)A_j(x, t)\big) + U(x, t)$$
We see that $H_C$ is already Weyl-ordered, so at least in this case we have $H = H_C$. On the other hand, the time evolution of the Wigner function can be written in an entirely tensorial way. This strongly suggests to me that $H = H_C$ is the natural choice for Euclidean space.
