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+d1)\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^{d1}$ when $k\rightarrow+\infty$ and thus the eigenvalue $k(k+d1)$ 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{d1}{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 LaplaceBeltrami 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}$?
Yes, the each gap survives a small perturbation (even a Hoelder perturbation), see here. For smoothness of the Laplacian with respect to the metric see here.
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.

$\begingroup$ 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. $\endgroup$ – Bazin Feb 2 '18 at 15:06
Indeed, each gap survives the small perturbation, but the smallness depends on the number of the gap. For generic metrics eigenvalues are distributed more uniformly.
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}$.
Actually, spectral gaps and clusters appear on manifolds with all geodesics closed.