Let $f \colon \mathbb{R}^3 \to \mathbb{R}^{2l+1} \,, l \in \mathbb{N}_0$ be an SO(3)-equivariant function, i.e.
$$f(\rho^1(g)r) = \rho^l(g)f(r) \,,\quad \forall g \in \text{SO(3)} \,, \forall r \in \mathbb{R}^3 \,,$$
where $(\rho^l, \mathbb{R}^{2l+1})$ is an irreducible representation of SO(3). Consider the subgroup $\text{S}(r) \subset \text{SO(3)}$ of rotations around the axis in direction of $r$ that leaves $r$ unchanged,
$$f(\rho^1(h)r) = f(r) = \rho^l(h)f(r) \,,\quad \forall h \in \text{S}(r) \,, \forall r \in \mathbb{R}^3 \,. \tag{1}\label{eq1}$$
We see that $f(r)$ is invariant under such a rotation. For $l=0$, $\rho^0$ is the trivial represenation and (1) holds automatically. For $l=1$, both $f(r)$ and $r$ are invariant under $\rho^1(h)$, thus they are parallel vectors $\in \mathbb{R}^3$. Thus, we see that $f$ is a radial function.
In general, (1) implies that $f(r)$ is an eigenvector of the map $\rho^l(h)$ with eigenvalue 1. I lack intuition about what this implies for the properties of $f$ in analogy to the $l=1$ case. I suspect there is a connection to the spherical harmonics, but I can't put my finger on it.
Can we draw conclusions or interpretations about the form of $f$ from (1) for $l>1$? Is there a connection to spherical harmonics?