Let $M$ be a compact manifold and let $U\subset M$ be an open set with a coordinate chart $\Phi:U\to V$ with $V\subset \mathbb{R}^d$. Suppose that $w\in C_c^\infty(U)$ is a smooth function supported in the open set $U$, and $F\in C^\infty(\mathbb{R})$ is another smooth function.
There are several ways one might sensibly quantize the semiclassical symbol $a_h(x,p):= w(x)F(|hp|^2)$. The first is to use a standard semiclassical quantization such as $$ {\operatorname*{Op}}_h(a_h) g(x) = (2\pi h)^{-d/2}\int e^{i\langle x-y,p\rangle/h} a_h(\Phi^{-1}(x), (D\Phi^{-1}_x)^Tp)g(y) dydp $$ for functions $g\in C^\infty(V)$ (or whatever is sensible, perhaps I made a mistake and the above is nonsensical).
The second way one might quantize $a_h$ is by setting directly $$ {\operatorname*{Op}}_h'(a) = w(x) F(h\Delta), $$ meaning for example that for an eigenfunction $\phi_j$ of the Laplacian, $\Delta\phi_j=\lambda_j\phi_j$, $$ ({\operatorname*{Op}}_h'(a)\phi_j)(y) := w(y) F(h\lambda_j)\phi_j(y). $$
My question is: to what extent should one expect ${\operatorname*{Op}}_h(a_h)$ and ${\operatorname*{Op}}'_h(a_h)$ to agree? To give a precise question: what assumptions on $F$ are necessary to ensure that these operators converge to each other in operator norm as $h\to 0$?
In general, what I am asking about is how one should go about comparing the spectral definition to the local definition. Any advice would be greatly appreciated.
There is a beautiful answer described by Michael Taylor in https://mtaylor.web.unc.edu/wp-content/uploads/sites/16915/2018/04/FIHA.pdf . The result that answers my question is Theorem 1.2, and for posterity I will sketch the idea below.
The first observation is that the hard part is just computing the principal symbol of $F(h^2\Delta)$, and in this case it is more convenient to consider $F(h^2\Delta)=\Phi(h\sqrt{-\Delta})$. Then one writes
$$ \Phi(h\sqrt{-\Delta}) = \int_{-\infty}^\infty e^{ish\sqrt{-\Delta}} \widehat{\Phi}(s) \,ds. $$
The operator $e^{ish\sqrt{-\Delta}}$ is the generator of the wave equation, and assuming $\Phi$ is smooth one only needs to consider $s\lesssim 1$ (because the Fourier transform should be rapidly decaying). Therefore the question of computing $\widehat{\Phi}(h\sqrt{-\Delta})$ reduces to integrals involving the wave equation for short time, the solution of which can be approximated using a Fourier integral operator.