How does one define the fractional Laplacian operator $(-\Delta)^s$ on a compact Riemannian manifold?
In $\mathbb{R}^n$, it is defined $$ (-\Delta)^s f(x) = c_{n,s} \int_{\mathbb{R}^n} \frac{f(x) - f(y)} {|x-y|^{n+2s}} \mathrm d y $$ where $c_{n,s}$ is a constant. Can we simply replace the integral to be an intergral over the manifold $M$ and change the denominator to $g(x,y)^{n+2s}$ where $g$ is the distance function?
No, that formula won't work. But you can use a spectral definition. There is an orthonormal basis $\{\phi_i\}_{i=1}^\infty$ for $L^2(M)$ consisting of eigenfunctions of $-\Delta$, with corresponding eigenvalues $\{\lambda_i\}$. Then for $u = \sum_i u_i \phi_i$, put $$(-\Delta)^s u(x) = \sum_{i=1}^\infty (\lambda_i)^s u_i \phi_i(x).$$