For example the eigenfunctions of a spherically symmetric membrane can be found in https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane. Then again I sometimes see people in Physics saying that the system is spherically symmetric and thus the solution also has to be spherically symmetric, for example here: https://youtu.be/DIuoFAW9H3E?feature=shared&t=583.
What is going on here? Under which circumstances can we say that the eigenfunctions of the Laplacian operator with spherically symmetric boundary conditions are themselves spherically symmetric?
The individual eigenfunctions have no reason to be invariant by rotation, only the eigenspaces. This is the case for your problem.
Indeed, say you have a symmetry $S$, and an operator $A$. The operator is invariant by the symmetry if it commutes with composition. Formally, for any $f$, $$ A(f\circ S) = (Af)\circ S $$ This is the case for the Laplacian with spherically symmetric boundary conditions and rotational symmetry. This only means that an eigenvector is mapped to an eigenvector with the same eigenvalue. If: $$ Af = \lambda f $$ then $$ A(f\circ S) = Af\circ S = \lambda f\circ S $$ so $f\circ S$ is still in the $\lambda$ eigenspace $E_\lambda$, therefore it is stable by composition via $S$. However, unless $\dim E_\lambda=1$, an $\lambda$ eigenvector has no reason to be invariant by the symmetry.
For example, in the wikipedia article, the eigenspace of $\lambda_{mn}$ is twice degenerate for $m\neq0$ due to the freedom to choose $C,D$. Therefore choosing an arbitrary basis for each of these 2D eigenspace, the rotation of a basis vector will still be a superposition of the original basis vectors.
Btw, I only considered scalar functions and symmetries being compostions, but in general, the functions can take values in an arbitrary vector space, and the symmetry can also act on the left as well.
Hope this helps.