Consider the case of Laplace's equation $$\nabla^2 f = 0$$ After the standard procedure of separation of variables, the general solution to this equation is found to be a linear combination of spherical harmonics $Y_\ell^m(\theta,\varphi)$: $$ f(r,\theta,\varphi) = \sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell f_\ell^m r^\ell Y_\ell^m (\theta,\varphi )$$
Surprisingly, the functions $Y_\ell^m(\theta,\varphi)$ also happen to be eigenfunctions of the operators $L_z$ and $\mathbf{L}^2$ where $$\mathbf{L}^2 = L_x^2+L_y^2+L_z^2$$ and the $L_x$, $L_y$ and $L_z$ operators are generators of the rotation group: $$ \begin{aligned} L_{x} &= \left(\sin\phi\frac{\partial}{\partial\theta}+\cot\theta\cos\phi\frac{\partial}{\partial\phi}\right) \\ L_{y} &= \left(-\cos\phi\frac{\partial}{\partial\theta}+\cot\theta\sin\phi\frac{\partial}{\partial\phi}\right) \\ L_{z} &= -\frac{\partial}{\partial\phi} \end{aligned} $$
This cannot be a coincidence as the $\nabla^2$ operator is invariant under rotations. This is something that I've noticed in many other PDEs solved by the separation of variables. Is there any deeper reason that explains why this happens?