Brief Background If you are familiar with Brownian motion in the sphere via SDEs in local coordinates/heat kernel, you can easily skip to the question. Otherwise, here is a brief background, which I hope is correctly written up.
Define Brownian motion on the sphere as the Markov process with transition density $p$ solving $$\frac{\partial p }{\partial t} = \frac12 \Delta_M p,$$ where $\Delta_M$ is the Laplace-Beltrami operator on the sphere. Using the parameterization, $$(\theta, \phi) \mapsto (\cos \phi \sin \theta, \sin \phi \sin \theta, \cos\theta)^T$$ and writing $\mathscr{L}=\frac12 \Delta_M$, we can write down the infinitesimal generator $\mathscr{L}$ in local coordinates $(\theta, \phi)$ as $$\mathscr{L}f = \frac12 \cot \theta \frac{\partial f}{\partial \theta} +\frac12 \left( \frac{\partial^2 f}{\partial \theta^2}+\frac{1}{\sin^2 \theta} \frac{\partial^2 f}{\partial \phi^2} \right).$$
If a steady-state $p_\infty(\theta, \phi)=\lim_{t \to \infty}p(t,\theta, \phi)$ exists, it must then be a harmonic function, i.e. $$\Delta_M p_\infty =0.$$
It is well known that the steady state PDF of the Langevin SDE $$dX_t = (\operatorname{div}_M g^{-1}-g^{-1} \nabla U)dt+\sqrt{2g^{-1}}dB_t$$ is $$p_\infty \propto e^{-U}\sqrt{\det g}$$ where $g$ is the metric tensor of the manifold, and $U$ is some potential defined in the local coordinates. Setting $U=0$, then the steady state of $$dX_t = (\operatorname{div}_M g^{-1})dt+\sqrt{2g^{-1}}dB_t$$ is $$p_\infty \propto \sqrt{\det g}.$$ Now this SDE only differs from the SDE for Brownian motion on a manifold by a factor of $1/2$ in the drift and $\sqrt{2}$ in the diffusion coefficient.
My Question For the sphere, we have $$g = \begin{pmatrix} 1 & 0 \\ 0 & \sin^2 \theta\\ \end{pmatrix},$$ so the steady state is therefore, $$p_\infty(\theta, \phi)= \sin(\theta),$$ since $\theta \in (0,\pi)$ for the parameterization above. I have verified indeed that this makes the RHS of the Fokker Planck equation associated to the last displayed SDE vanish. However, $$\Delta_M p_\infty \neq 0,$$ indeed we obtain $$\Delta_M \sin \theta = \frac{\cos^2 \theta-\sin^2\theta}{\sin\theta} \neq 0.$$ What can explain this discrepancy? In other words, I have a derived a steady state PDF but it is not a harmonic function on the manifold, which appears to be a contradiction. Am I missing something simple, or have I made a mistake in understanding the theory?
Thanks in advance and let me know if I should include anymore details in the background or the question.