I have always found this kind of things very confusing myself, but now that I have some idea I'll try to address it. You can find a proof of a stronger statement in Proposition 7.7 in [this][1] work of mine with Luca Gennaioli, but part of the proof is hidden in Proposition 7.5. Note that we only assume that $(M,g)$ is stochastically complete, and this is not even really needed (see the **Note** below). The main point is the following. Let me call
$$ (-\Delta)^{s}_{\rm B} \, u := \frac{1}{\Gamma(-s)} \int_{0}^{\infty} (e^{t\Delta}u-u)\frac{dt}{t^{1+s}}$$
the Bochner's fractional Laplacian, this may seem unrelated but for me it is the key to understand how this works. Moreover, note that for $\lambda \ge 0$ and $s \in (0,1)$ there holds 
$$\varphi(\lambda) := \lambda^{s} =  \frac{1}{\Gamma(-s)} \int_0^\infty
(e^{-\lambda t}-1)\frac{dt}{t^{1+s}} . $$
This is an equality between real numbers and has nothing to do with the Laplacian. Then, from here spectral theory does most of the work. 

I claim the following fact, which answers your question in great generality.

**Claim.** Let $X \subset L^2(M)$ be a linear subspace with the following property: whenever $u \in X$ then $(-\Delta)_B^s \, u $ exists in the sense of Bochner (that is, the right-hand side in its definition converges absolutely in $L^2(M)$). Then $X\subseteq {\rm Dom}((-\Delta)^{s}_{\rm Spec})$ and
$$ (-\Delta)_{\rm Spec}^s u := \int_{\sigma(-\Delta)} \lambda^{2s} d\langle E_\lambda u, \cdot \rangle = (-\Delta)_B^s \, u \,.$$

Note that the last equality alone implies $X\subseteq {\rm Dom}((-\Delta)^{s}_{\rm Spec})$ as
$$\int_{\sigma(-\Delta)} \lambda^{2s} d\langle E_\lambda u, u \rangle = \| (-\Delta)_{\rm Spec}^s u  \|^2_{L^2(M)} =  \| (-\Delta)_B^s \, u \|^2_{L^2(M)} \\= \left\| \int_0^\infty (e^{t\Delta}u-u) \frac{dt}{t^{1+s/2}} \right\|_{L^2(M)}^2 \le  \left( \int_0^\infty \| e^{t\Delta}u-u \|_{L^2(M)} \frac{dt}{t^{1+s/2}} \right)^2 <+\infty$$
since $u\in X$ and we have used Minkowski's Integral Inequality. 

Basically you need the property in the hypothesis in the claim with $X=C_c^\infty(M)$ and in my work with Luca we prove it for $X=H^{2s+\epsilon}(M)$ (for every $\epsilon \ll 1$) on every stochastically complete $(M,g)$. I believe this to be true with $H^{2s}(M)$ as it is know in bounded domains $\Omega \subset \mathbb{R}^n$, but in a setting so general as every complete Riemannian manifold one needs to be careful. Every proof I have seen of this for $\Omega \subset \mathbb{R}^n$ heavily uses the discreteness of the spectrum and
interpolation theory, and Fourier transform the proofs for the entire $\mathbb{R}^n$. 

The fact that the hypothesis is true for $X=C_c^\infty(M)$ is easy, you just need to prove that the integral in $(-\Delta)_B^s \, u$ is absolutely convergent near $0$ and $\infty$. At the origin just using $\|e^{t\Delta} u-u \|_{L^2(M)} \le t \|u\|_{C^2}(M)$ gives convergence, and at $\infty$ use that $ \|e^{t\Delta} u-u \|_{L^2(M)} \le \|e^{t\Delta} u \|_{L^2(M)} + \|u \|_{L^2(M)} \le 2\|u \|_{L^2(M)} $. 

It only remains to prove the claim, and the proof goes as follows. For every $v\in L^2(M)$ by standard spectral theory (see e.g. equation (A.49) in A. Grigor'yan's book "Heat Kernel and Analysis on Manifolds" for why is the first equality true) 
$$ \langle \varphi(-\Delta) u, v \rangle_{L^2} = \int_{0}^\infty \varphi(\lambda) d\langle E_\lambda u,v \rangle = \frac{1}{\Gamma(-s)}  \int_0^\infty \left( \int_0^\infty
(e^{-\lambda t}-1)\frac{dt}{t^{1+s}} \right) d\langle E_\lambda u,v \rangle \,.$$
Since the integrand is always negative by Tonelli's theorem you can exchange the order of integration and you get
$$ \langle \varphi(-\Delta) u, v \rangle_{L^2} = \int_0^\infty \langle e^{t\Delta}u-u,v \rangle_{L^2} \frac{dt}{t^{1+s}} = \langle (-\Delta)_B\, u,v \rangle_{L^2} <+\infty \,, $$
and this is finite and all the integrals converge since $u\in X$ satisfies the hypothesis in the claim.


**Note:** indeed we assume $(M,g)$ being stochastically complete because we want, at the same time, to prove that the Bochner's fractional laplacian and the Singular Integral fractional Laplacian coincide, but you are uninterested in this. The fact that Bochner's fractional Laplacian is in $L^2(M)$ is still true even without stochastical completeness, and this what you really need.

  [1]: https://arxiv.org/pdf/2306.11590.pdf