The (real) spherical harmonics $Y^m_\ell$, where where $\ell=0,1,2,\dots$ and $m = -\ell, -\ell+1,\dots,\ell-1,\ell$, are a set of eigenfunctions of the Laplace-Betrami operator $\Delta_{\mathbb S^2}$ on the sphere that satisfy
\begin{align}
\Delta_{\mathbb S^2} Y^m_\ell = -\ell(\ell+1)Y^m_\ell. \tag{1}
\end{align}
They form an orthonormal basis of the Hilbert space $L^2(\mathbb S^2)$ of square-integrable (real-valued) functions, which means that we can expand any $f\in L^2(\mathbb S^2)$ as
\begin{align}
f = \sum_{\ell=0}^\infty\sum_{m=-\ell}^\ell f_{\ell m} Y^m_\ell, \tag{2}
\end{align}
with $f_{\ell m}\in\mathbb R$. Now consider $\Delta_{\mathbb S^2} f$.
Question: under what conditions on $f$ can we interchange the infinite sum and the Laplace-Beltrami operator? In other words, when do we have the following?
\begin{align}
\Delta_{\mathbb S^2}f \stackrel{?}{=} \sum_{\ell=0}^\infty\sum_{m=-\ell}^\ell f_{\ell m} \Delta_{\mathbb S^2} Y^m_\ell = -\sum_{\ell=0}^\infty\sum_{m=-\ell}^\ell f_{\ell m} \ell(\ell+1)Y^m_\ell \tag{3}
\end{align}
In particular, is it sufficient to require that $f$ be infinitely differentiable? If so, this would be a quick way of showing that there exist no (smooth) harmonic functions on the sphere other than the constant functions.
EDIT: And furthermore, is there a way to deduce from the coefficients $f_{\ell m}$ wether $f$ is sufficiently differentiable to even apply the Laplace-Beltrami operator to it?
EDIT: The answer by Giuseppe Negro is great, but I haven't been able to convince myself why the first theorem in the answer should be true, and I also could not find a reference proving it. I think I just found an alternative proof, though, which seems even simpler. I would love to hear if people agree with it. Here it is.
Claim: (3) holds if $f$ is smooth (i.e. infinitely differentiable).
To simplify notation, condense the subscripts $\ell,m$ on $Y^\ell_m$ into a single subscript $i$, denote the corresponding eignenvalues by $\lambda_i$, and denote the $L^2$ inner product between two functions $f,g$ by $\langle f,g\rangle$.
Proof. Let $f\in C^\infty(\mathbb S^2)\subset L^2(\mathbb S^2)$ and write $f = \sum a_iY_i$ with $a_i\in\mathbb R$.Then we also have $\Delta_{\mathbb S^2}f\in C^\infty(\mathbb S^2)\subset L^2(\mathbb S^2)$ and hence $\Delta_{\mathbb S^2}f$ can be written as $\Delta_{\mathbb S^2}f = \sum b_iY_i$ with $b_i\in\mathbb R$, where coefficient $b_j$ is given by $b_j=\langle\Delta_{\mathbb S^2}f, Y_j\rangle$. The crucial property of $\Delta_{\mathbb S^2}$ is that is symmetric with respect to its action on smooth functions (this follows from Green's second identity), i.e. since both $f$ and each $Y_j$ are smooth, we have
\begin{align}
b_j&=\langle\Delta_{\mathbb S^2}f, Y_j\rangle = \langle f, \Delta_{\mathbb S^2} Y_j\rangle = \langle f, \lambda_j Y_j\rangle \\
%
&= \lambda_j \langle f, Y_j\rangle = \lambda_j a_j.
\end{align}
This shows that $\Delta_{\mathbb S^2}f = \Delta_{\mathbb S^2}(\sum a_iY_i) = \sum \lambda_i a_i Y_i = \sum a_i \Delta_{\mathbb S^2} Y_i$, as desired.
Does this work or am I overlooking something?
Best Answer
The general theorem is very simple and not specific of the Laplace--Beltrami. You can even formulate a version of the following in abstract functional analytic terms, for closed self-adjoint operators in Hilbert space.
GENERAL OPERATOR-THEORETIC FRAMEWORK.
For a general $f\in L^2(\mathbb S^{d-1})$ write $\widehat{f}(\ell, m)$ to denote the coefficient in $$ f=\sum_{\ell\ge 0}\sum_m \widehat{f}(\ell, m) Y_{\ell, m}, $$ where $Y_{\ell, m}$ denotes a complete orthonormal system of spherical harmonics of degree $\ell$. Note that requiring $f\in L^2(\mathbb S^{d-1})$ is exactly the same as requiring $\sum_{\ell, m}\lvert \widehat{f}(\ell, m)\rvert^2<\infty$. Finally, recall that $-\Delta_{\mathbb S^{d-1}}Y_{\ell, m}=\ell(\ell+d-2)Y_{\ell, m}$.
Theorem. Let $f\in L^2(\mathbb S^{d-1})$ be such that $\sum_{\ell\ge 0} \ell(\ell+d-2)\sum_m\lvert \widehat{f}(\ell, m)\rvert^2<\infty$. Then $f\in \mathrm{Dom}(-\Delta_{\mathbb S^{d-1}})$ and $$\tag{1} -\Delta_{\mathbb S^{d-1}} f = \sum_{\ell\ge 0}\sum_m \ell(\ell+d-2)\widehat{f}(\ell, m) Y_{\ell, m}, $$ where the series converges in $L^2(\mathbb S^{d-1})$ sense.
This can actually be taken as the definition of $\mathrm{Dom}(-\Delta_{\mathbb S^{d-1}})$, in the sense of closed self-adjoint operators.
REFORMULATION OF THE QUESTION IN SOBOLEV SPACE TERMINOLOGY.
Your question now reads as follows. Is it true that $C^\infty(\mathbb S^{d-1})\subset \mathrm{Dom}(-\Delta_{\mathbb S^{d-1}})$? You can also ask: in this case, is the series in (1) uniformly convergent?
I bet the answer is affirmative to both questions, motivated by the analogy with usual Fourier series. Of course the usual method of proof for Fourier series, which is based on integration by parts, is not immediately available here. (But see the end of this post for a spherical version of that argument).
REFERENCE. See Section 3.8: "Sobolev spaces on the sphere" in the book of Atkinson and Han: https://link.springer.com/book/10.1007/978-3-642-25983-8
Note that Atkinson and Han define norms $$ \lVert f\rVert_{H^s(\mathbb S^{d-1})}^2:=\sum_{\ell\ge 0}(\ell+\delta_d)^{2s}\sum_m \lvert \widehat{f}(\ell, m)\rvert^2,$$ where $\delta_d=\frac{d-2}{2}$, but this exact value is not important for us here. They say that $f\in H^s(\mathbb S^{d-1})$ if and only if $\lVert f\rVert_{H^s(\mathbb S^{d-1})}<\infty$. With this terminology, $$ f\in \mathrm{Dom}(-\Delta_{\mathbb S^{d-1}}) \iff f\in H^2(\mathbb S^{d-1}).$$ Now as it is briefly explained in the book, $C^\infty(\mathbb S^{d-1})\subset H^s(\mathbb S^{d-1})$ for all $s\ge 0$ (and actually it is dense). In particular, $C^\infty(\mathbb S^{d-1})\subset H^2(\mathbb S^{d-1})$. $\Box$
DIRECT ANSWER AND CONCLUSION. It seems that the Atkinson and Han reference is also not satisfactory. So let us prove directly that $f\in C^\infty(\mathbb S^{d-1})$ implies that $\lvert\widehat{f}(\ell, m)\rvert\le C_N (1+\ell)^N$ for each $N\in\mathbb N$, where $C_N>0$ may depend on $f$ but it does not depend on $\ell$. Once this is proven, the fact that $C^\infty(\mathbb S^{d-1})\subset H^2(\mathbb S^{d-1})$ is immediate.
The key of this proof is integration by parts. By definition, since $Y_{\ell, m}$ is a orthonormal system, $$ \widehat{f}(\ell, m)=\int_{\mathbb S^{d-1}} f(y)\overline{Y_{\ell, m}(y)}\, d\sigma(y).$$ Now, $$Y_{\ell, m}=\frac{(-\Delta_{\mathbb S^{d-1}})^k Y_{\ell, m}}{(\ell(\ell+d-2))^k}, $$ so \begin{equation*} \begin{split} \lvert \widehat{f}(\ell, m)\rvert &\le\frac{1}{(\ell(\ell+d-2))^k} \left\lvert \int_{\mathbb S^{d-1}} (-\Delta_{\mathbb S^{d-1}})^k f(x) \overline{Y_{\ell, m}(x)}\, d\sigma(x)\right\rvert\\ &\le \frac{\lVert (-\Delta_{\mathbb S^{d-1}})^k f\rVert_{L^2(\mathbb S^{d-1})}}{(\ell(\ell+d-2))^k} \end{split} \end{equation*}
and now the claim follows immediately, since $f\in C^\infty(\mathbb S^{d-1})$ implies that the norm in the right--hand side of the latter inequality is finite for all $k$. $\Box$
OUTDATED REFERENCE. In a previous version of this answer I cited the book of Feng Dai and Yuan Xu: https://link.springer.com/book/10.1007/978-1-4614-6660-4