Do Laplace-Beltrami eigenfunctions vary continuously with the metric? I'm interested in Laplace Beltrami operators $$-\Delta_g:\ \   D(-\Delta_g) \longrightarrow L^2\left(M,\sqrt{|g|}dx\right)$$
on a smooth compact Riemannian Manifold (M,g). Let us fix a unique metric $g$ on $M$. 
For any other smooth metric $\widetilde g$ on $M$, we can identify the square integrable functions with respect to its associated volume form with our original $L^2$ above, via the unitary map
$U:f\longmapsto \sqrt{|g|/|\widetilde g|}f.$
Under this unitary identification the Laplace-Beltrami operators corresponding to the various smooth metrics on $M$ yield a family of operators on a common domain of definition.
It has been proved that the eigenvalues of this family depend continuously on the metrics.
Is this true for the eigenfunctions or eigenprojections as well?
I feel like the answer should be yes, arguing in coordinates that if the coefficients of the operators are close, their resolvents and hence their spectral projections are, but I am not sure how to prove it rigorously.
I know that if I consider an analytic one parameter family of metrics the result holds essentially by Kato's perturbation theory.
But I am ultimately interested in proving that a certain composition of maps from functions on the manifold to $\mathbb R$  is robust to arbitrary, small (in the $C^\infty$-topology) changes in the metric and one of the maps in the composition is an operator of the form $f(-\Delta_g)$.
This is my first question here and also a duplicate of this question on stackexchange, so if I'm transgressing against any etiquette rules, please do let me know. :)
 A: EDITED: Added clarification, as pointed out by @TerryTao.
Let $g_1$, $g_2$ be Riemannian metrics and $\Delta_1$, $\Delta_2$ their respective Laplacians. Let $\lambda_1$ be an eigenvalue of $\Delta_1$ and $\lambda_2$ an eigenvalue of $\Delta_2$.
I think what can be proved is the following: Given an eigenfunction $u_2$ of $\Delta_2$ with eigenvalue $\lambda_2$, there exists an eigenfunction $u_1$ of $\Delta_1$ with eigenvalue $\lambda_1$ such that
$$ \|u_2-u_1\|_2 \le C(|\lambda_2-\lambda_1| + \|g_2-g_1\|_{2,2})\|u_2\|_2, $$
where $C$ depends on both $g_1$ and $\lambda_1$ (specifically, the spectral gap between $\lambda_1$ and the other eigenvalues of $\Delta_1$), $\|\cdot\|_{2,2}$ is the $W^{2,2}$ Sobolev norm with respect to a background metric $g_0$, and $\|\cdot\|_2$ is the $L^2$ norm with respect to $dV_1$.
Here's a sketch of my proof:
Let $\lambda_2$ be an eigenvalue of $\Delta_2$ and $u_2$ a nontrival eigenfunction. Let $\lambda_1$ be an eigenvalue of $\Delta_1$ and $u_1$ a nontrival eigenfunction.
A straighforward calculation shows that
$$
  (-\Delta_1+\lambda_1)(u_2-u_1)
  = (\Delta_2-\Delta_1)u_2 - (\lambda_2-\lambda_1)u_2.
$$
Choose $u_1$ so that, for any eigenfunction $u$ of $\Delta_1$ with eigenvalue $\lambda_1$,
$$ \int u(u_2-u_1)\,dV_1 = 0. $$
The claim now follows by elliptic estimates.
