Gaps in the spectrum of Laplace-Beltrami operators

Let us consider $\mathbb S^d$ the unit Euclidean sphere of $\mathbb R^{d+1}$ and let $\Delta_{\mathbb S^d}$ be the Laplace operator on $\mathbb S^d$. We have $$-\Delta_{\mathbb S^d}=\sum_{k\in \mathbb N}k(k+d-1)\mathbb P_k, \quad I=\sum_{k\in \mathbb N}\mathbb P_k,$$ where $\mathbb P_k$ is the orthogonal projection on $\mathcal S_k$, the spherical harmonics with degree $k$. The dimension of $\mathcal S_k$ is equivalent to $c_d k^{d-1}$ when $k\rightarrow+\infty$ and thus the eigenvalue $k(k+d-1)$ has a large multiplicity when $d\ge 2$. Looking at $\sqrt{-\Delta_{\mathbb S^d}}$, we find that the eigenvalues are $$\{0\}\cup_{k\ge 1}\{\lambda_k=k\bigl(1+\frac{d-1}{k}\bigr)^{1/2}\},$$ corresponding to the eigenspace $\mathcal S_k$. As a result, we have gaps in the spectrum of $\sqrt{-\Delta_{\mathbb S^d}}$ since $$\lambda_{k+1}-\lambda_k\sim 1, \quad k\rightarrow+\infty.$$ Now my question: if we perturb smoothly the metric on $\mathbb S^d$, we get a new Laplace-Beltrami operator $\tilde{\Delta}$. Are the gaps in the spectrum of $\sqrt{-\Delta_{\mathbb S^d}}$ surviving to the perturbation ? Are they still present for $\sqrt{-\tilde\Delta}$?

But maybe, Weyl's asymptotic formula (see p155 of Chavel: Eigenvalues in Riemannian Geometry) $$(\lambda_k)^{d/2} \sim \frac{(2\pi)^d k}{\text{Vol}(D^d).\text{Vol}(M)}$$ (which holds for each compact Riemannian manifold of dimension $d$) is sufficient for you.
• Thanks a lot, I am indeed very grateful for your help. I do not believe that Weyl's law is enough for the stability of gaps since it probably could be the same with eigenvalues spread on the interval $(k,k+1)$ instead of being clustered at the integers. Also, it is surprising to me that the references that you kindly indicate are so recent; I should say that I thought that, for a smooth perturbation it was either wrong or classical. Commented Feb 2, 2018 at 15:06
On the other hand, if Laplacian on the sphere is perturbed by a lower order operator, eigenvalues will be still divided into eigenvalue clusters separated by spectral gaps. The asymptotic width of such clusters depends on the order of the perturbation and other things. Say, for $$-\Delta +V(x)$$ the width of the cluster is $$\lesssim \lambda^{-1/2}$$, but it is $$\lesssim \lambda ^{-3/2}$$ if $$V(x)=-V(-x)$$, while the width of the spectral gap is $$\asymp \lambda^{1/2}$$.