The question says it all. I can elaborate in the morning.
This is ultimately related to question Can a well-behaved, positive-definite function $\phi(x)$ always be represented by $\phi(x)=\sum_{n=0}^{\infty}\mid a_{n}\psi(x)_n\mid^{2} $?
I plan on using both notions in concert.
spherical harmonics are eigenfunctions of the angular part of the laplace operator. and they happen to be complete and ortonormal. so a function of three variables can always be expanded as: $$ f(x,y,z) \overset{\longrightarrow}{\textrm{sc}} f(r,\theta,\phi)=\sum_{l,m}f_{l,m}(r)Y_l^m(\theta,\phi) $$ where "sc" means a change from $(x,y,z)$ to spherical coordinates $(r,\theta,\phi)$
in 3 dimensions we have: $$ \Delta=\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}-\frac{\Lambda^2}{r^2} $$ where $\Lambda^2$ is a differential operator involving only angles. now if we take a polynomial of G-degree with the following properties: $$ p_G(\lambda\vec{x})=\lambda^G\, p_G(\vec{x}) $$ $$ \Delta p_G =0 $$ then one could define: $$ Y_G(\vec{\omega})=r^{-G}p_G(\vec{x})\,\,\,\, ,\,\,\,\,|\vec{x}|=r \,\,\, , \,\,\, |\vec{\omega}|=1 $$ if you calculate: $$ 0=\Delta p_G=\Delta {r^G Y_G}=...=r^{G-2}\big(G(G+1)-\Lambda^2 \big)Y_G $$ we find: $$ \Lambda^2\, Y_G=G(G+1) Y_G $$ so this demonstrates that an homogeneous harmonic polynomial is an eigenfunction of the angular part of the laplacian, by placing further conditions one can actually obtain the expression of a spherical harmonic.