I saw in many physics texts term by term differentiation of spherical harmonics expansion, but since they're physics texts they're without rigourous proof. Take for example the following from Wikipedia
Given the multipole expansion of a scalar field
${\displaystyle \phi =\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\phi _{\ell m}(r)Y_{\ell m}(\theta ,\phi ),}$
we can express its gradient in terms of the VSH as
${\displaystyle \nabla \phi =\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\left({\frac {d\phi _{\ell m}}{dr}}\mathbf {Y} _{\ell m}+{\frac {\phi _{\ell m}}{r}}\mathbf {\Psi } _{\ell m}\right).}$
I know the spherical harmonics $Y_{\ell m}$ forms an orthonormal basis for $L^2$ functions, so that any $L^2$ function $\phi$ can be written as an infinite series expansion of spherical harmonics in the above way. If $\phi$ is $H^1$ function, then we expect its derivative to also have an spherical harmonic expansion. However, how do we know rigourosuly that it can be obtained by term by term differentiation of the series for $\phi$ in the above manner? In other words, how do we know the coefficients for the spherical harmonic expansion of $\nabla\phi$ is related to the coefficients for the spherical harmonic expansion of $\phi$ via term by term differentiation? Since the convergence of the series that I know of is only in $L^2$, not uniform, I can't quite see how to properly justify it.. Could one justify it if we have more regularity like $C^k$ for some large $k$?
The sad truth is that this will not be true in general. However, there is a magic wand, called distribution theory. Your series converges in the distributional sense and there we have the freshman‘s dream come true— convergence is preserved under differentiation. This is shown in any standard text (fortunately, there are elementary approaches which avoid the use of functional analysis, in particular, duality theory for locally convex spaces).