Using the classic spherical harmonics theory, one obtains the $k$-th eigenvalue of the $n$-dimensional round sphere $S^n$ to be $k(k+n-1)$, and its multiplicity is $\binom{n+k}{k}-\binom{n+k-1}{k-1}$, see e.g. [Berger, Gauduchon,Mazet, "Le spectre d'une variété riemannienne", Lecture Notes in Mathematics, Vol. 194 Springer-Verlag]. By looking at eigenfunctions of the Laplacian on $S^n$,$S^{2n+1}$ and $S^{4n+3}$ (note they are the unit spheres of $\mathbb R^{n+1}$, $\mathbb C^{n+1}$ and $\mathbb H^{n+1}$) that are respectively invariant under the natural actions of $\mathbb Z_2$, $S^1$ and $S^3$, one can obtain the eigenfunctions hence the $k$-th eigenvalue of the projective spaces $\mathbb R P^n$, $\mathbb C P^n$ and $\mathbb H P^n$ respectively. These are, respectively, $2k(n+2k-1)$, $4k(n+k)$ and $4k(k+2n+1)$.

QUESTION: Compute the multiplicity of the $k$-th eigenvalue $4k(k+2n+1)$ of the Laplacian on $\mathbb H P^n$.

The multiplicity of the $k$-th eigenvalue on the sphere arises as the dimension of the space of homogeneous harmonic polynomials of degree $k$ in $\mathbb R^{n+1}$, which is the difference between the dimensions of the spaces of homogeneous polynomials of degree $k$ and $k-1$ in $\mathbb R^{n+1}$. This follows from a polar decomposition of the first space as direct sum of spaces of homogeneous polynomials. The multiplicity of the $k$-th eigenvalue $4k(n+k)$ of the Laplacian on the complex projective space $\mathbb C P^n$ is given by the difference of the squares of these dimensions, it is $\binom{n+k}{k}^2-\binom{n+k-1}{k-1}^2$. This again follows from a polar decomposition of the space of harmonic homogeneous polynomials on $\mathbb C^{n+1}$. However, I do not know how to extend this idea to decompose the space of harmonic homogeneous polynomials on $\mathbb H^{n+1}$ and hence obtain the multiplicity of the $k$-th eigenvalue of $\mathbb H P^n$. This actually seems to be an exercise of [Berger, Gauduchon, Mazet]. Any natural guesses with similar differences of binomial powers seem to be wrong, since the multiplicity of the first eigenvalue (i.e., $k=1$) of $\mathbb H P^n$ should be $2n^2+3n$.

HINT: Apparently the desired multiplicity should coincide with the dimension of the vector space formed by homogeneous polynomials $f$ in $\mathbb H^{n+1}=\mathbb R^{4n+4}$ with coordinates $(x_\alpha,y_\alpha,z_\alpha,w_\alpha), 1\leq\alpha\leq n+1$, of degree $2k$ such that $Lf=0$, where $L=\sum_j L_jL_j$ and
$$L_1=\sum_\alpha y_\alpha\frac{\partial}{\partial x_\alpha}-x_\alpha\frac{\partial}{\partial y_\alpha}+w_\alpha\frac{\partial}{\partial z_\alpha}-z_\alpha\frac{\partial}{\partial w_\alpha}
$$
$$L_2=\sum_\alpha z_\alpha\frac{\partial}{\partial x_\alpha}-w_\alpha\frac{\partial}{\partial y_\alpha}-x_\alpha\frac{\partial}{\partial z_\alpha}+y_\alpha\frac{\partial}{\partial w_\alpha}
$$
$$L_3=\sum_\alpha w_\alpha\frac{\partial}{\partial x_\alpha}+z_\alpha\frac{\partial}{\partial y_\alpha}-y_\alpha\frac{\partial}{\partial z_\alpha}-x_\alpha\frac{\partial}{\partial w_\alpha}
$$
The space of homogeneous polynomials of degree $2k$ in $4n+4$ variables has dimension $\binom{4n+2k+3}{2k}$, so the desired multiplicity should be this number minus the dimension of the subspace formed by polynomials such that $Lf=0$. Note that the above operators $L_j$, when restricted to the unit sphere $S^{4n+3}$ of $\mathbb H^{n+1}$ give the three vector fields that span the distribution tangent to the Hopf fibers $S^3$.