The shell theorem and the Hairy Ball theorem To get straight to the question, skip to the 'This is what spherical symmetry really means' paragraph.
I recently inquired about the shell theorem but this is a different matter (link), so I hope it is ok to make a new post. There is an elegant way of deriving the shell theorem by taking advantage of the spherical symmetry of the situation. But instead of just saying "by symmetry", I was pondering how exactly to describe that precisely. I have some thoughts but would like to know what more experienced people think.
So suppose $\rho: \mathbb{R}^3 \to \mathbb{R}$ is some spherically symmetric mass distribution; spherically symmetric means for any $T \in \operatorname{SO}(3)$ and $r \in \mathbb{R}^3$ one has $\rho(Tr) = \rho(r)$. The gravitational field at a point $r$ is
$$ g(r) = -G \int \rho(s) \frac{r - s}{|r - s|^3} dV $$
so that at the point $Tr$ we have
$$ g(Tr) = -G \int \rho(s) \frac{Tr - s}{|Tr - s|^3} dV. $$
Now we can apply the change of variables theorem. Since $\det T = 1$ it gives us
$$ g(Tr) = -G \int \rho(Ts) \frac{Tr - Ts}{|Tr - Ts|^3} dV$$
Note $\rho(Ts) = \rho(s)$, and $|Tr - Ts| = |r - s|$ since $T$ is a rotation, so we can pull $T$ out of the integral to get
$$g(Tr) = T \left( - G\int \rho(s) \frac{r-s}{|r-s|^3}dV \right) = Tg(r).$$
This is what spherical symmetry really means. Now let $\Sigma$ be some sphere centered around the origin, and let's say large enough that all of the mass $\rho$ lies inside of it. Then knowning that $gT = Tg$ does immediately give us that the magnitude of $g$ is the same everywhere in $\Sigma$, since rotations preserve magnitude. But it does not seem to immediately tell us, from a mathematical standpoint at least, that $g$ points inward everywhere.
I thought of a way out: in $\Sigma$ you can decompose the field $g = g_{\_} + g_{\perp}$ into its tangent and normal components (and the action of $\operatorname{SO}(3)$ respects this decomposition). Now $g_{\_}$ is a tangent vector field on the 2-sphere $\Sigma$. By the Hairy ball theorem it has a zero somewhere. But since rotations carry 0 to 0, it must be 0 everywhere!  So $g$ only has a perpendicular component.
This is nice, but it only works because we are in an odd-dimensional space. So would the shell theorem not work for gravity in a plane? Is there a way to argue this not involving the hairy ball theorem?
 A: Theorem. Let $n > 2$ be a positive integer. Assume that in the space
$\mathbb{R}^n$ (or possibly in $\mathbb{R}^n \setminus \{0\}$) we
have a vector field $g(r)$ with the property that for any $T \in
\text{SO}(n)$
$$T\,g(r) = g(Tr)$$
Then $g(r)$ is radial, i.e. there exists a scalar function $w :
\mathbb{R}_{+} \to \mathbb{R}$ such that
$$g(r) = w(|r|)\, r$$
Furthermore, if we define the function $$V(r) = W(|r|) \,\, \text{ where } \,\, W(\lambda) = \int_{\lambda_0}^{\lambda} \,s\, w(s) \, ds$$
then for any $r\in \mathbb{R}^n \setminus \{0\} $
$$\nabla \, V(r) = g(r)$$ i.e. all rotationally-invariant vector fields on $\mathbb{R}^n \setminus \{0\}$ are both radial and potential fields.
$$ $$
Proof: Let $r \in \mathbb{R}$ be an arbitrary non-zero vector. Define the satbilizer of $r$ from the rotation group $\text{SO}(n)$
$$\text{Stab}(r) = \{\, T \in \,\text{SO}(n) \, : \, T\,r = r\}$$ Furthermore, denote by $L$ the orthogonal complement of the vector $r$, which is by definition
$$L = \{\, v \in  \mathbb{R}^n \, : \, v \cdot r = 0 \,\}$$
Then $\dim L = n-1$ and it is $\text{Stab}(r)$ invariant, i.e. for any $v \in \mathbb{R}^n$ such that $v \cdot r = 0$ and for any $T \in \text{Stab}(r)$
\begin{align*}
T\,v \cdot r = T\,v \cdot T\, r = v \cdot r = 0
\end{align*}
so $T\,v \in L$. Decompose the space $$\mathbb{R}^n = L \oplus \mathbb{R} \,r$$
Then $v \in \mathbb{R}^n$ decomposes uniquely as $v = v_L + \lambda \,r$, with $v_L \in L$.
For any $T \in \text{Stab}(r)$
$$T\, v = T(v_L + \lambda \,r) = T\, v_L + \lambda \,T \, r = w_L + \lambda \, r$$ where $w_L = T\,v_L \in L$.
Since $\text{Stab}(r)$ is a subgorup of $\text{SO}(n)$ and as such preserves the dot product on $\mathbb{R}^n$, when restricted to $L$ it preserves
the dot product on the $n-1$ dimensional subspace $L$. Moreover, if you take any linear transformation $T_L$ on $L$ that preserves the dot product,
it can be extended to $T \in \text{SO}(n)$ as $$T \,v = T\, (v_L + \lambda\, r) = T_L\, v_L +  \lambda\, r$$ Thus, $\text{Stab}(r)$ is the full
rotation group of the $n-1$ dimensional subspace $L$ and is therefore isomorphic to $\text{SO}(n-1)$. Form the properties of $\text{SO}(n-1)$ we know
that if $v \in L$ is any non-zero vector, then its orbit under $\text{Stab}(r)$ is an $n-2$ dimensional sphere, so
there is always a transformation $T \in \text{Stab}(r)$ such that $T \,v \neq v$.
Now, let us focus on the $\text{SO}(n)$-invariant vector field $g(r)$. Take any non-zero vector $r \in \mathbb{R}^n$. Split the space
$$\mathbb{R}^n = L \oplus \mathbb{R} \,r$$ as explained above. Then $g(r)$ decomposes uniquely into
$$g(r) = g_L + w\, r$$
where $g_L \in L$ and $w \in \mathbb{R}$. For any $T \in \text{Stab}(r)$
$$T\,g(r) = g(T\,r) = g(r)$$ which in decomposed form translates into
$$T\,g(r) = T(g_L + w\, r) = T\,g_L + w\, T\,r = T\,g_L + w\, r = g_L + w\,r$$ and when we cancel out the term $w\,r$ from both sides of the last identity,
we find out that for any $T \in \text{Stab}(r)$
$$T\,g_L = g_L$$
But this is possible if and only if $g_L = 0$, because as I already mentioned before,
for a non-zero vector $g_L$ from $L$, there is always a transformation $T \in \text{Stab}(r)$ such that $T \,g_L \neq g_L$
(because the orbit is a proper $n-2$ dimensional sphere).
Hence, for any non-zero $r \in \mathbb{R}^n$ there exists a real number $u(r) \in \mathbb{R}$, that varies with respect to $r$,
such that $$g(r) = u(r)\, r$$ i.e. the vector field is radial.
Now, by invariance, for any $T \in \text{SO}(n)$ and for any non-zero $r \in \mathbb{R}^n$,
$$T (u(r)\, r) = u(r)\, T\,r = T\,g(r) = g(T\,r) = u(T\,r) \, T\,r$$ which, because of the identity $u(r)\, T\,r = =  u(T\,r) \, T\,r$,
is possible if and only if $u(r) = u(T\,r)$.
Fix one unit vector $r_0 \in \mathbb{R}^n$ ($|r_0| = 1$). Take any other non-zero vector $r \in \mathbb{R}^n$. Both vectors $\frac{r}{|r|}$
and $r_0$ lie on the unit sphere and since $\text{SO}(n)$ acts transitively on the unit sphere (in fact with stabilizers always isomorphic to $\text{SO}(n-1)\,$)
there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $$T \, \frac{r}{|r|} = r_0$$ Now do the following:
\begin{align}
u(r) =& u(T\,r) \\
u(r)\, =& \, u\left(\,T \left(|r| \, \frac{r}{|r|}\right)\,\right) \\
u(r) \, =&\, u\left(|r| \, T \, \frac{r}{|r|}\right) \\
u(r) \, =&\, u\big(\,|r| \, r_0\,\big)
\end{align}
It follows from the last identity that there exists a single-variable function $$w : \mathbb{R}_{+} \to \mathbb{R}$$
$$w(\lambda) = u(\lambda \, r_0)$$
such that
$$u(r) = w(|r|)$$
All of this allows us to conclude that for $n > 2$, any $\text{SO}(n)$-invariant vector field $g(r)$ on $\mathbb{R}^n$ is radial and can be expressed as
$$g(r) = w(|r|)\, r $$ for some suitable function $w : \mathbb{R}_{+} \to \mathbb{R}$. If we set
$$W(\lambda) = \int_{\lambda_0}^{\lambda} \, s\, w(s)\,ds$$
$$V(r) = W(|r|)$$ then $\frac{d}{d\lambda} W(\lambda) = W'(\lambda) = \lambda\, w(\lambda)$ and hence by the chain rule
$$\nabla\, V(r) = \nabla \, W(|r|) = W'(|r|) \, \frac{r}{|r|} = |r|\, w(|r|)  \, \frac{r}{|r|} =  w(|r|)\, r = g(r)$$ i.e.
the rotationally-invariant vector field $g(r)$ is potential.
$$ $$
$$ $$
Older Version.
Assume that for an odd dimensional space $\mathbb{R}^n$ we have a vector field $g(r)$ with the property that for any $T \in \text{SO}(n)$
$$T\,g(r) = g(Tr)$$
Theorem. Let $n$ be odd. For every non-zero vector $r \in \mathbb{R}^n$ there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $Tr = r$ and the space of all vectors such that $Tv = v$ is one dimensional, spanned by $r$. In other words, for any vector $r \in \mathbb{R}^n$ there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $r$ is an egienvector of $T$ with eigenvalue $1$, where the multiplicity of the eigenvalue is $1$.
Proof: Take an arbitrary non-zero vector $r \in \mathbb{R}^n$ and define the unit vector $\hat{r} = \frac{r}{|r|}$. Denote by $L$ the $n-1$ dimensional vector subspace of $\mathbb{R}^n$ orthogonal to $\hat{r}$. Since $L$ is non-trivial, choose a unit vector $\hat{r}_2 \in L$. Then $\hat{r}$ and $\hat{r}_2$ are pairwise orthogonal. Continue this procedure. Assume you have selected $k$ unit vectors $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k$, that are pairwise orthogonal. Then, there exists a unique $n-k$ dimensional vector subspace $L_k$ which is the orthogonal complement of the span of $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k$. If $k < n$ then $L_k$ is non-trivial and you can find aunit vector $\hat{r}_{k+1} \in L_k$. Thus, you have found $k+1$ unit vectors $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k, \, \hat{r}_{k+1}$, that are pairwise orthogonal. Continue  like this until $k = n$. You have constructed an orthonormal basis $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_n$ of $\mathbb{R}^n$. Think of this new vectors as column-vectors. Define the orthogonal matrix
$$U = \big[ \hat{r} \, \, \hat{r}_2 \,\, ... \, \hat{r}_n \big] \, \in \, \text{SO}(n)$$ If you denote by $e_1 = [1, \, 0, \, ..., \, 0]^T$  the first basis vector of $\mathbb{R}^n$, interpreted as a column-vector, then
$$\hat{r} = U \, e_1 \,\, \text{ and therefore } \,\, e_1 = U^T \, \hat{r}$$
If $n = 2\,m + 1$, construct the block-diagonal matrix
$$R = R(\theta_1,..., \theta_m) = \begin{bmatrix} 1 &  & & &\\
                                & R_2(\theta_1) & & & &\\
& & R_2(\theta_2)  & & \\& &  & ... & \\ & &  & & R_2(\theta_m) 
\end{bmatrix}$$
where for each $j=1,..,m$
$$R_2(\theta_j) = \begin{bmatrix} \cos(\theta_j) &   -\,\sin(\theta_j)\\
                                \sin(\theta_j)  &   \cos(\theta_j)
\end{bmatrix}$$ and $\theta_j \in (0, 2\pi)$. By construction, $$R(\theta_1,..., \theta_m)\, v = v \,\, \text{ if and only if } \,\, v = \lambda \,e_1$$, i.e. this matrix has eigenvalue $1$ with multiplicity $1$. The rest of the eigenvalues are true complex unitary numbers. Finally, define the orthogonal transformation
$$T = U \, R \, U^T \in \text{SO}(n)$$
By construction, $$T \, r = T \, |r| \, \frac{r}{|r|} = |r|\, T \,\hat{r} = 
 |r|\, U\, R \, U^T \,\hat{r} = |r|\, U\, R\, e_1 =  |r|\, U\, e_1 = |r|\, \hat{r} = r$$ Moreover, since $T$ and $R$ are conjugate matrices, they have the same eigenvalues, including multiplicities. Therefore $T$ has eigenvalue $1$ of multiplicity $1$ and eigenvector $r$. End of proof
Now, you can apply this theorem to any $\text{SO}(n)$-invariant vector field on $\mathbb{R}^n \setminus \{0\}$, i.e. $$g : \mathbb{R}^n \setminus \{0\} \to \mathbb{R}^n$$
$$g(T\,r) = Tg(r) \,\,\text{ for any } \,\, T \in \text{SO}(n)$$ For an arbitrary non-zero vector $r \in \mathbb{R}^n$, pick one orthogonal transfromation $T_r \in \text{SO}(n)$ with the property described in the theorem above. Hence $T_r\, r = r$. Furthermore
$$T_rg(r) = g(T_rr) = g(r)$$ which means that $g(r)$ is also an eigenvector of $T_r$ with egienvalue $1$. Since the eigenspace of $T_r$ that corresponds to the eigenvaule $1$ is one dimensional and spanned by $r$, the vector $g(r)$ must be a multiple of $r$, i.e. there exists a scalar $\lambda(r) \in \mathbb{R}$ such that $$g(r) = \lambda(r) \, r$$ As you can see, the vector field $g(r) = \lambda(r) \, r$ is radial.
What happens if you take $n=2$? Can you find a counterexample? (I think it is fairly straightforward)
Now, there is another approach to your question, in particular that if you define the potential function
$$V(r) = G \, \int_{\Sigma} \, \frac{\rho(s)}{|r - s|} \, ds$$ and the function $\rho(s)$ is nice enough so that we can differentiate with respect to $r$ under the integral, then
$$g(r) = \nabla \,V(r) = - \, G \, \int_{\Sigma} \, \rho(s)\,\frac{(r-s)\,\,}{|r - s|^3} \, ds$$
So basically, your $\text{SO}(n)$-invariant vector field is in fact a potential field and comes from a potential function $V(r)$ which is also $\text{SO}(n)$-invariant, i.e. for any $T \in \text{SO}(n)$
$$V(Tr) = V(r)$$ Fix one unit vector $r_0 \in \mathbb{R}^n$ ($|r_0| = 1$). Now take any other non-zero vector $r \in \mathbb{R}^n$. Both vectors $\frac{r}{|r|}$
and $r_0$ lie on the unit sphere and since $\text{SO}(n)$ acts transitively on the unit sphere (in fact with stabilizers always isomorphic to $\text{SO}(n-1)\,$) there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $$T \, \frac{r}{|r|} = r_0$$ Now do the following:
\begin{align}
V(r) =& V(T\,r) \\
V(r)\, =& \, V\left(\,T \left(|r| \, \frac{r}{|r|}\right)\,\right) \\
V(r) \, =&\, V\left(|r| \, T \, \frac{r}{|r|}\right) \\
V(r) \, =&\, V\big(\,|r| \, r_0\,\big)
\end{align}
It follows from the last identity that there exists a single-variable function $$W : \mathbb{R}_{+} \to \mathbb{R}$$
$$W(\lambda) = V(\lambda \, r_0)$$
such that
$$V(r) = W(|r|)$$ Then, by the chain rule, the gradient of $V$ with respect to $r$ is
$$\nabla \, V(r) = \nabla \, W(|r|) =  \frac{dW}{d\lambda}(|r|)\, \frac{r}{|r|} = W'(|r|)\, \frac{r}{|r|}$$
i.e. the vector field $g(r) = \nabla \, V(r) = W'(|r|)\, \frac{r}{|r|}$ is radial.
