If you are only interested in expanding polynomials then you can forgo Hilbert spaces and get an exact expansion into a finite linear combination.
Spherically-invariant case Let $f$ be a polynomial in $x_1,\ldots, x_n$ (the coefficients may be real, complex, or even a field of characteristic zero). If $f$ is invariant under isometries of the $n$-dimensional space, i.e. $f(Ax)=f(x)$ for any $n\times n$ orthogonal matrix $A,$ then $f=g(r^2),$ where $g$ is a polynomial in one variable of degree $\deg(f)/2$ and
$$r^2=x_1^2+x_2^2+\ldots+x_n^2.$$
This is intuitively clear: a rotationally invariant function depends only on the distance to the center and $g(t)=f(t,0,\ldots,0)$ is a polynomial. (For $n\geq 2$ it's sufficient that $f$ be invariant under rotations, i.e. the special orthogonal group $SO(n)$.)
The theory of spherical harmonics, which is closely related to the representation theory of the orthogonal group $O(n),$ generalizes this observations to polynomials (and even more general functions) that are not necessarily spherically invariant.
General case The first key theorem states that a polynomial $f$ can be expanded as
$$ f(x)=\sum_{k=0}^d r^{2k}h_k(x) \qquad (*)$$
where $d\leq\deg(f)/2$ and $h_k$ is a harmonic polynomial of degree at most $\deg(f)-2k.$ Recall that a polynomial is $h$ is harmonic if $\Delta h=0,$ where $\Delta$ is the Laplace operator,
$$\Delta h=\frac{\partial^2 h}{\partial x_1^2}+\frac{\partial^2 h}{\partial x_2^2}+\ldots+\frac{\partial^2 h}{\partial x_n^2}.$$
Moreover, the harmonic polynomials $h_k$ in $(*)$ are uniquely determined. If $f$ is spherically invariant then $h_k$ is simply the coefficient of $t^k$ in the polynomial $g(t)$ from above, viewed as a constant (i.e. degree $0$) polynomial.
So now every polynomial can be expanded in terms of the powers of $r^2$ and harmonic polynomials. What can we say about the latter? Let $\mathcal{H_\ell}$ (respectively, $\mathcal{H}$) denote the space of homogeneous harmonic polynomials of degree $\ell$ (respectively, all harmonic polynomials). Every harmonic polynomial is uniquely decomposed into the sum of homogeneous polynomials of various degrees, and these components are themselves harmonic:
$$\mathcal{H}=\bigoplus_{\ell\geq 0}\mathcal{H}_\ell\qquad (**)$$
The second key theorem states that $\mathcal{H}_\ell$ is an irreducible representation of the orhtogonal group $O(n)$ of dimension $\binom{n+\ell-1}{n-1}-\binom{n+\ell-3}{n-1}$ (for $\ell=0$ or $1$, the second term is zero). For $n\geq 3,$ this representations remains irreducible upon restriction to the special orthogonal group $SO(n)$ (for $n=3$, this is the rotation group).
For small $n$, namely $n=2,3,4$, there are standard bases in the spaces $\mathcal{H}_\ell$ that have been tabulated and extensively studied. In particular, for $n=3,$ $\mathcal{H}_\ell$ is a $(2\ell+1)$-dimensional representation of $SO(3)$, or spin $\ell$ representation in physics language, which has a basis indexed by an integer $m,\ -\ell\leq m\leq \ell$ constructed using associate Legendre polynomials $P_\ell^m.$ Depending on your goals, you may want a suitable explicit description of these functions in rectangular or spherical coordinates (the Wikipedia article Spherical harmonics, including the references, is a good starting point).
Let me also mention a common point of confusion related to your question 3: in view of $(*),$ for any polynomial $f$ as above, there exists a unique harmonic polynomial $h$ such that the restrictions of $f$ and $h$ to the unit sphere
$$S^{n-1}: x_1^2+x_2^2+\ldots+x_n^2=1$$ coincide, $f|_{S^{n-1}}=h|_{S^{n-1}}.$ In other words, the space $\mathcal{H}$ of harmonic polynomials can be naturally identified with the space of polynomial functions on the unit sphere $S^{n-1}$ and the decomposition $(**)$ becomes the decomposition of the polynomial functions on the sphere into irreducible representations of $O(n),$ each of which occurs with multiplicity one. These functions are frequently called spherical harmonics.
Addendum Here is a high level view of polynomial spherical harmonics from a representation theory vantage point. The decomposition $(*)$ is related to the representation theory of the Lie algebra $\mathfrak{sl}_2.$ The operators
$$E=r^2/2,\ F=-\Delta/2,\ H=\deg+n/2$$
on the vector space of polynomials in $n$ variables commute with orthogonal transformations and form a representation of $\mathfrak{sl}_2.$ (I've described the skew-symmetric analogue in this answer.) Homogeneous harmonic polynomials are precisely lowest weight vectors for $\mathfrak{sl}_2,$ and the second key theorem amounts to saying that the lowest weight spaces are irreducible representations of $O(n).$ One consequence of this description is that the coefficients $h_k(x)$ in $(*)$ can be found inductively starting with $h_d(x)$ using repeated applications of the Laplace operator (the precise statement is omitted due to bulky notation). Furthermore, the first and second key theorems can be combined into the statement that the space $\mathcal{P}$ of polynomials in $n$ variables has the following decomposition into irreducible components under the joint $O(n)$ and $\mathfrak{sl}_2$ actions:
$$\mathcal{P}=\bigoplus_{\ell\geq 0} \mathcal{H}_\ell \otimes V_\ell, $$
where $V_\ell$ is the lowest weight $\mathfrak{sl}_2$-module with lowest weight $\ell+n/2.$
This is one of the starting points of Roger Howe's theory of reductive dual pairs. In the present case, the reductive dual pair is $(O(n),SL(2,\mathbb{R}))$ over the real numbers.
Remark The group $SL(2)$ is secretely a symplectic group, in fact, the isometry group of a 2-dimensional symplectic vector space, and the theory extends to reductive dual pairs consisting of an orthogonal group and a symplectic group of any rank. The same theory over $\mathbb{Q}$ and its adeles is behind some classical results in the theory of theta functions due to Siegel and Weil. In particular, it explains the importance of theta functions with harmonic coefficients.
