Let $Y_n^j, \, -n\leq j \leq n$ be a spherical harmonic. Let $f$ be an even function on the sphere. Then, $f$ can be introduced as $$ f=\sum_{j=-n}^A\hat f(j)Y_n^j, $$ where $\sum_{j=-n}^A|a_j|>0$ and $A$ is the largest integer in $[-n, n]$ such that $\hat f(A)\neq 0$.
Now, consider spherical polynomial (finite combination of spherical harmonics). We can just take degree of the function $f$, say $f^2$ $$ f^2=\sum_{j=-n}^A\hat f(j)^2(Y_n^j)^2+2\sum_{-n\leq j<l\leq A}\hat f(j)\hat f(l)Y_n^jY_n^l $$ Question: spherical polynomial is the finite combination of spherical harmonics. Can one take limit there?
Yes, as you anticipate, although the spirit of "every function can be representated as an infinite sum of spherical harmonics", there are problems with convergence, as well as problems with the meaning of multiplication of two infinite sums of functions. These issues already arise (famously) on the one-sphere and the theory of Fourier series.
As a clear simple case: if the function $f$ is smooth, then the infinite sum of spherical harmonics converges pointwise uniformly to it, can be differentiated termwise, and all derivatives' infinite series expansions converge uniformly pointwise as well. And then it is immediate that the series can be squared and express $f^2$.
Similarly, if $f$ is merely sufficiently differentiable, its infinite expression in spherical harmonics converges (absolutely) uniformly pointwise, and can be squared as you indicated.
Details about how much differentiability, of various sorts, suffices are more complicated.
EDIT: For references for rigorous treatment: Stein and Weiss' old book "Introduction Fourier Analysis on Euclidean Spaces" has a chapter that covers many set-up issues about rigorous treatment of spherical harmonics, including things like estimate of sup norms in terms of $L^2$ norms. Also, Wikipedia's page on "spherical harmonics" does include several remarks about convergence arguments. My own course notes from introductory modular forms courses (at http://www.math.umn.edu/~garrett/m/mfms/) include some essays that treat the $L^2$ Sobolev theory, in the context of Weyl's criterion for equidistribution on spheres: http://www.math.umn.edu/~garrett/m/mfms/notes_2014-15/09_spheres.pdf
One basic assertion about convergence via Sobolev theory is the following. For a test function $f$ on the sphere, and non-negative integer $k$, $|f|^2_k=\langle (1-\Delta)^k f,f\rangle$ is the $k$th Sobolev norm (squared). Let $H^k$ be the completion of test functions with respect to that norm. Then a basic (not-so-trivial) theorem is that for $k>n/2$ the space $H^k$ lies inside the space of continuous functions, and the convergence of the finite partial sums of the spherical harmonic expansion converge uniformly pointwise to the function.