I) Let us work in units where $G=1=\rho$.
In this answer, we will additionally assume that the massive object is star-shaped. We can then use spherical coordinates $(r,\theta,\varphi)$. The surface profile is then given as
$$\tag{1} {\bf r}~=~r{\bf n}~=~f({\bf n}){\bf n}, \qquad r~=~f({\bf n})~=~\sum_{\ell m} c_{\ell m} Y_{\ell m}({\bf n})~\geq~0 ,\qquad c_{\ell m}~=~\int_{S^2}\!d^2n ~Y^{\ast}_{\ell m}({\bf n})f({\bf n}),$$$$ {\bf r}~=~r{\bf n}~=~f({\bf n}){\bf n}, \qquad r~=~f({\bf n})~=~\sum_{\ell m} c_{\ell m} Y_{\ell m}({\bf n})~\geq~0 ,$$ $$ c_{\ell m}~=~\int_{S^2}\!d^2n ~Y^{\ast}_{\ell m}({\bf n})f({\bf n}),\tag{1}$$
where ${\bf n}\in S^2$ is a unit vector. In eq. (1) we have used DumpsterDoofus' idea to expand in spherical harmonics. The volume
$$\tag{2} V[f]~:=~ \int_V \! d^3r ~=~\frac{1}{3}\int_{S^2}\!d^2n ~ f({\bf n})^3 .$$$$ V[f]~:=~ \int_V \! d^3r ~=~\frac{1}{3}\int_{S^2}\!d^2n ~ f({\bf n})^3 .\tag{2}$$
is a functional of the surface profile $f({\bf n})$.
$$\tag{3} \delta V ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}),\qquad \delta^2 V ~=~2\int_{S^2}\!d^2n ~ f({\bf n})\delta f({\bf n})^2.$$$$ \delta V ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}),\qquad \delta^2 V ~=~2\int_{S^2}\!d^2n ~ f({\bf n})\delta f({\bf n})^2.\tag{3}$$
The volume is kept fixed by a constraint
$$\tag{4} V[f]~=~V_0.$$$$ V[f]~=~V_0.\tag{4}$$
Potential:
$$\tag{5} -\Phi({\bf r})~:=~\int_{V} \! \frac{d^3r^{\prime}}{|{\bf r}-{\bf r}^{\prime}|} ~=~\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\! \frac{r^{\prime 2}dr^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|} .$$$$ -\Phi({\bf r})~:=~\int_{V} \! \frac{d^3r^{\prime}}{|{\bf r}-{\bf r}^{\prime}|} ~=~\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\! \frac{r^{\prime 2}dr^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|} .\tag{5}$$
$$\tag{6} -\delta \Phi({\bf r}) ~=~\int_{S^2}\!d^2n^{\prime} \frac{f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|{\bf r}-f({\bf n}^{\prime}){\bf n}^{\prime}|}.$$$$ -\delta \Phi({\bf r}) ~=~\int_{S^2}\!d^2n^{\prime} \frac{f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|{\bf r}-f({\bf n}^{\prime}){\bf n}^{\prime}|}.\tag{6}$$
Field:
$$\tag{7} -{\bf g}({\bf r})~:=~ \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{{\bf r}-{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3} .$$$$ -{\bf g}({\bf r})~:=~ \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{{\bf r}-{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3} .\tag{7}$$
Radial field:
$$ -g_r({\bf r})~=~\frac{\partial\Phi({\bf r})}{\partial r} ~=~{\bf n}\cdot \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{r-{\bf n}\cdot{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3}$$ $$~=~\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\!r^{\prime 2}dr^{\prime} \frac{r-{\bf n}\cdot{\bf n}^{\prime}r^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|^3}.\tag{8}$$$$ \begin{align}-g_r({\bf r})~=~&\frac{\partial\Phi({\bf r})}{\partial r} ~=~{\bf n}\cdot \nabla \Phi({\bf r}) ~=~\int_{V} \!d^3r^{\prime} \frac{r-{\bf n}\cdot{\bf r}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|^3}\cr ~=~&\int_{S^2}\!d^2n^{\prime}\int_0^{f({\bf n}^{\prime})}\!r^{\prime 2}dr^{\prime} \frac{r-{\bf n}\cdot{\bf n}^{\prime}r^{\prime}}{|{\bf r}-r^{\prime}{\bf n}^{\prime}|^3}.\end{align}\tag{8}$$
II) Potential energy:
$$ U[f]~:=~ -\iint_{V\times V} \! \frac{d^3r~ d^3r^{\prime}}{2|{\bf r}-{\bf r}^{\prime}|} ~=~ \frac{1}{2}\int_V \! d^3r ~\Phi({\bf r}) $$ $$~=~ -\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime} \int_0^{f({\bf n})}\int_0^{f({\bf n}^{\prime})} \frac{r^2dr~ r^{\prime 2}dr^{\prime}}{2|r{\bf n}-r^{\prime}{\bf n}^{\prime}|}.\tag{9}$$$$ \begin{align} U[f]~:=~& -\iint_{V\times V} \! \frac{d^3r~ d^3r^{\prime}}{2|{\bf r}-{\bf r}^{\prime}|} ~=~ \frac{1}{2}\int_V \! d^3r ~\Phi({\bf r}) \cr ~=~& -\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime} \int_0^{f({\bf n})}\int_0^{f({\bf n}^{\prime})} \frac{r^2dr~ r^{\prime 2}dr^{\prime}}{2|r{\bf n}-r^{\prime}{\bf n}^{\prime}|}.\end{align}\tag{9}$$
$$\tag{10} \delta U ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n})~\Phi(f({\bf n}){\bf n}).$$$$ \delta U ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n})~\Phi(f({\bf n}){\bf n}).\tag{10}$$
$$ \delta^2 U~=~-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}$$ $$ -\int_{S^2}\!d^2n~ f({\bf n})^2 ~\delta f({\bf n})^2g_r(f({\bf n}){\bf n}) +2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \Phi(f({\bf n}){\bf n}).\tag{11}$$$$ \begin{align}\delta^2 U~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}\cr &-\int_{S^2}\!d^2n~ f({\bf n})^2 ~\delta f({\bf n})^2g_r(f({\bf n}){\bf n})\cr &+2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \Phi(f({\bf n}){\bf n}).\end{align}\tag{11}$$
III) To treat the constraint (4) we use Lagrange multiplier method. We should minimize the functional
$$\tag{12} E[f]~:=~U[f] +\lambda (V_0-V[f]).$$$$ E[f]~:=~U[f] +\lambda (V_0-V[f]).\tag{12}$$
The first variation is $$\tag{13} \delta E ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}) \{\Phi(f({\bf n}){\bf n})-\lambda\}.$$$$ \delta E ~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n}) \{\Phi(f({\bf n}){\bf n})-\lambda\}.\tag{13}$$
We conclude that
The surface potential $\Phi(f({\bf n}){\bf n})=\lambda<0$ of a stationary shape $f$ is a constant, i.e. independent of the unit vector ${\bf n}$.
The second variation is
$$ \delta^2 E~=~-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}$$ $$ -\int_{S^2}\!d^2n~ f({\bf n})^2 \delta f({\bf n})^2g_r(f({\bf n}){\bf n}) +2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \underbrace{\{\Phi(f({\bf n}){\bf n})-\lambda\}}_{=0}.\tag{14}$$$$ \begin{align}\delta^2 E~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~ \frac{f({\bf n})^2\delta f({\bf n})~ f({\bf n}^{\prime})^2\delta f({\bf n}^{\prime})}{|f({\bf n}){\bf n}-f({\bf n}^{\prime}){\bf n}^{\prime}|}\cr &-\int_{S^2}\!d^2n~ f({\bf n})^2 \delta f({\bf n})^2g_r(f({\bf n}){\bf n})\cr &+2\int_{S^2}\!d^2n~f({\bf n})\delta f({\bf n})^2 \underbrace{\{\Phi(f({\bf n}){\bf n})-\lambda\}}_{=0}.\end{align}\tag{14}$$
For systematic reasons, the second variations (14) should be consistent with the constraint (4) to first order.
IV) We finally consider a ball $f({\bf n})=R$. We then have
$$V~=~\frac{4\pi}{3}R^3, \qquad {\bf g}({\bf r})~=~-\frac{4\pi}{3}{\bf r}, \qquad \Phi({\bf r})~=~2\pi\left(\frac{r^2}{3}-R^2\right), \qquad \Phi(R)~=~-\frac{4\pi}{3}R^2, \qquad U~=~-\frac{16\pi^2}{15}R^5.\tag{15}$$$$V~=~\frac{4\pi}{3}R^3, \qquad {\bf g}({\bf r})~=~-\frac{4\pi}{3}{\bf r}, \qquad \Phi({\bf r})~=~2\pi\left(\frac{r^2}{3}-R^2\right), $$ $$ \Phi(R)~=~-\frac{4\pi}{3}R^2, \qquad U~=~-\frac{16\pi^2}{15}R^5.\tag{15}$$
We want to show that the ball is a stable stationary shape. We have to show that the second variation (14) is semipositive definite. Eq. (3) simplifies to
$$\tag{16} \frac{\delta V}{R^3} ~=~\int_{S^2}\!d^2n ~\delta f({\bf n}) ~=~\sqrt{4\pi}\delta c_{00},$$$$ \frac{\delta V}{R^3} ~=~\int_{S^2}\!d^2n ~\delta f({\bf n}) ~=~\sqrt{4\pi}\delta c_{00},\tag{16}$$
$$\tag{17}\frac{\delta^2 V}{2R^3}~=~\int_{S^2}\!d^2n ~|\delta f({\bf n})|^2 ~=~\sum_{\ell m}|\delta c_{\ell m}|^2.$$$$\frac{\delta^2 V}{2R^3}~=~\int_{S^2}\!d^2n ~|\delta f({\bf n})|^2 ~=~\sum_{\ell m}|\delta c_{\ell m}|^2.\tag{17}$$
Since the constraint (4) should be maintained to first (but not necessarily second) order, we must demand that the zeromode vanishes
$$\tag{18}\delta c_{00}~=~0.$$$$\delta c_{00}~=~0.\tag{18}$$
The second variation (14) simplifies to (see e.g. Ref. 1)
$$ \frac{\delta^2 E}{R^3} ~=~-\int_{S^2}\!d^2n \int_{S^2}\!d^2n^{\prime}~\frac{\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}}{|{\bf n}-{\bf n}^{\prime}|} -\frac{g_r(R)}{R}\int_{S^2}\!d^2n~\delta f({\bf n})^2 $$ $$~=~-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}\sum_{\ell m} \frac{4\pi}{2\ell+1}Y^{\ast}_{\ell m}({\bf n})Y_{\ell m}({\bf n}^{\prime})$$ $$+\frac{4\pi}{3}\int_{S^2}\!d^2n~|\delta f({\bf n})|^2 $$ $$~=~4\pi\sum_{\ell m}\left(\frac{1}{3}-\frac{1}{2\ell+1}\right)|\delta c_{\ell m}|^2 ~\geq~0,\tag{19}$$$$ \begin{align}\frac{\delta^2 E}{R^3} ~=~&-\int_{S^2}\!d^2n \int_{S^2}\!d^2n^{\prime}~\frac{\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}}{|{\bf n}-{\bf n}^{\prime}|} -\frac{g_r(R)}{R}\int_{S^2}\!d^2n~\delta f({\bf n})^2 \cr ~=~&-\int_{S^2}\!d^2n\int_{S^2}\!d^2n^{\prime}~\delta f({\bf n})~\delta f({\bf n}^{\prime})^{\ast}\sum_{\ell m} \frac{4\pi}{2\ell+1}Y^{\ast}_{\ell m}({\bf n})Y_{\ell m}({\bf n}^{\prime})\cr &+\frac{4\pi}{3}\int_{S^2}\!d^2n~|\delta f({\bf n})|^2 \cr ~=~&4\pi\sum_{\ell m}\left(\frac{1}{3}-\frac{1}{2\ell+1}\right)|\delta c_{\ell m}|^2 ~\geq~0,\end{align}\tag{19}$$
which is non-negative since the zeromode (18) is absent. So the ball is a stable stationary shape.
References:
- J.D. Jackson, Classical Electrodynamics; chapter 3.