Do Laplace-Beltrami eigenfunctions vary continuously with the metric?

Question

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. :)

Answer 1

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.