I) Let us work in units where $G=1=\rho$.

In this answer, we will additionally assume that the massive object is [star-shaped](http://en.wikipedia.org/wiki/Star_domain). We can then use [spherical coordinates](http://en.wikipedia.org/wiki/Spherical_coordinate_system) $(r,\theta,\varphi)$. The surface profile is then given as

$$ {\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](http://en.wikipedia.org/wiki/Spherical_harmonics). The volume 

$$ 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})$. 

$$ \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 

$$ V[f]~=~V_0.\tag{4}$$

Potential:

$$ -\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}$$

$$ -\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: 

$$ -{\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:

$$ \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:

$$ \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}$$

$$ \delta U
~=~\int_{S^2}\!d^2n ~ f({\bf n})^2\delta f({\bf n})~\Phi(f({\bf n}){\bf n}).\tag{10}$$

$$ \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](http://en.wikipedia.org/wiki/Lagrange_multiplier). We should minimize the functional

$$   E[f]~:=~U[f] +\lambda (V_0-V[f]).\tag{12}$$ 

The first variation is
$$ \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

$$ \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), $$
$$ \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 

$$ \frac{\delta V}{R^3} ~=~\int_{S^2}\!d^2n ~\delta f({\bf n})
~=~\sqrt{4\pi}\delta c_{00},\tag{16}$$ 

$$\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

$$\delta c_{00}~=~0.\tag{18}$$

The second variation (14) simplifies to (see e.g. Ref. 1)

$$ \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:

1. J.D. Jackson, _Classical Electrodynamics;_ chapter 3.