Let us set ourselves in $\mathbb R^d$ and let $\Omega$ be an open set. Let $P$ be a differential operator with coefficients that admits a fundamental solution $E$ which is $C^\infty$ outside of $0$. If $f \in C^\infty (\Omega)$ then we know that the solution $u$ of the equation $Pu=f$ is $C^\infty$ in $\Omega$.
My question is if we can do the same for the Laplace-Beltrami operator on the sphere (or even in a more general manifold) to show (without explicit calculations) that the solutions of $-\Delta + \lambda$ are $C^\infty$ on the sphere?
The reason for asking this is that in papers on nodal sets, the space that the equation is given on is not specified. I assume it is $L^2$ in order to apply the spectral theorem, but am wondering if the same phenomena of bounded open sets on $\mathbb R^d$ appear on the compact Riemannian manifolds too.
Yes, elliptic regularity still applies on manifolds. The $C^\infty$ smoothness is a local property, so it suffices to work within a coordinate patch. The change of variable by a coordinate map transforms a uniformly elliptic PDE with smooth coefficients into another uniformly elliptic PDE with smooth coefficients, but on a Euclidean domain. Hence the results. For references, see Elliptic Regularity on Manifolds.