The general formula to get the potential energy of any spherical distribution is this :
\begin{equation}\tag{1}
U = – \int_0^R \frac{GM(r)}{r} \, \rho(r) \, 4 \pi r^2 \, dr,
\end{equation}
where $M(r)$ is the mass inside a shell of radius $r < R$. It is easy to get the gravitational energy of a uniform sphere of mass $M$ and radius $R$ :
\begin{equation}\tag{2}
U = -\, \frac{3 G M^2}{5 R}.
\end{equation}
In general, for any spherical distribution of total mass $M$ and exterior radius $R$, we can write this :
\begin{equation}\tag{3}
U = -\, \frac{k \, G M^2}{R},
\end{equation}
where $k > 0$ is a constant that depends on the internal distribution. $k = \frac{3}{5}$ for the uniform distribution. For a thin spherical shell of radius $R$ (all mass concentrated on its surface), we can get $k = \frac{1}{2}$.
Now, I suspect that for all cases :
\begin{equation}\tag{4}
\frac{1}{2} \le k < \infty.
\end{equation}
Physically, this makes sense. But how to prove this from the general integral (1) ?
To simply things a bit, we may introduce the dimensionless variable $x = r/R \le 1$, and defines relative mass $\bar{M}(x) \equiv M(r)/M \le 1$ and relative density $\bar{\rho}(x) = \rho(r) / \rho_{\text{average}}$, where $\rho_{\text{average}} = 3 M/4 \pi R^3$. Thus, integral (1) takes the following form :
\begin{equation}\tag{5}
U = -\, \frac{3 G M^2}{R} \int_0^1 \bar{M}(x) \, \bar{\rho}(x) \, x \, dx.
\end{equation}
The last integral is $\frac{k}{3}$. I'm not sure this may help to prove (4).
Best Answer
Start by writing that: $$\tag1 M=\int^R_0 \rho(r)4\pi r^2 dr = \int^R_0 \frac{dM}{dr}dr.$$ From (1), we can identify that: $$\tag2 \frac{dM}{dr}dr = 4\pi r^2 \rho(r)dr.$$ Now let's return to your equation (1): $$U = - \int_0^R \frac{GM(r)}{r} \, \rho(r) \, 4 \pi r^2 dr,$$ using equation (2), this becomes: $$U=-G\int_0^R \frac{M(r)}{r} \frac{dM}{dr}dr.$$
Let's concentrate on just the integral part, and call it $I$ for concreteness, $$I=\int_0^R \tag{3}\frac{M(r)}{r} \frac{dM}{dr}dr.$$ Aha! This looks like something we can integrate by parts. Let's set $u=\frac{M(r)}{r}$, and $v'=\frac{dM}{dr}$. Then $$u'=-\frac{M(r)}{r^2} + \frac{1}{r}\frac{dM}{dr}.$$ Continuing with IBP, we get: $$I=\Bigg[\frac{M(r)}{r}M(r)\Bigg]^R_0-\int^R_0 \Big(-\frac{M(r)}{r^2} + \frac{1}{r}\frac{dM}{dr}\Big)M(r) dr.$$ Let's look at solely the integral part, calling it $J$: $$J=-\int^R_0 \frac{M^2(r)}{r^2} dr + \underbrace{\int^R_0 \frac{M(r)}{r}\frac{dM}{dr} dr}_{\textrm{we know this integral, it's }I}.$$ So now we have that: $$I=\Bigg[\frac{M^2(r)}{r}\Bigg]^R_0 + \int^R_0 \frac{M^2(r)}{r^2} dr - I,$$ and so: $$2I= \Bigg[\frac{M^2(r)}{r}\Bigg]^R_0 + \int^R_0 \Bigg(\frac{M(r)}{r}\Bigg)^2 dr.$$ For the case of non-pathological $M$, that is $M(0)=0$, we get: $$I=\frac{M^2}{2R} + \frac{1}{2}\int^R_0 \Bigg(\frac{M(r)}{r}\Bigg)^2 dr.$$
Hence we can identify from this that $U$ has the limit of $-\frac{GM^2}{2R}$, and from that we subtract the integral of something positive definite.
I believe this proves (to physicist standards) your conjecture about $k$.
Request of simplifying the process:
Let's just consider with no dimensions the integral: $$K=\int^X_0 \frac{f(x) f'(x)}{x} dx.$$
Take: $$g(x)=\frac{f^2(x)}{2x},$$ then taking the derivative: $$g'(x) = \frac{f(x) f'(x)}{x}-\frac{f^2(x)}{2 x^2}.$$ We can recognise our original integral in this mix, and so integrating $g'(x)$ will give us our integral, less $\int^X_0 \frac{f^2(x)}{2 x^2} dx$, and so we can write: $$K=\big[g(x)\big]^X_0 + \frac{1}{2}\int^X_0 \frac{f^2(x)}{x^2}dx= \Bigg[\frac{f^2(x)}{2 x}\Bigg]^X_0 + \frac{1}{2}\int^X_0 \bigg( \frac{f(x)}{x} \bigg)^2 dx.$$ Giving us the result we wanted. As IBP is always the inverse of the product rule, with a little bit of thought you can generally figure out the original function. However if I'd just posted this last 4-step proof, I think you'd have just thought I was a wizard...