What are the gravitational binding energies of giant planets?

Question

What are the gravitational binding energies of the planets in our solar system? In particular, interested in the giant planets: Jupiter, Saturn, Uranus, and Neptune. Ideally the information would be from a paper or other peer-reviewed source, but I'll take what I can get. What I'd like to avoid is the approximation $$ U_\circ=\frac{3G}{5}\cdot\frac{m^2}{r}=\frac{16\pi^2G}{15}r^5\rho^2 $$ which comes from assuming that the planet is a sphere of uniform density. This isn't all that realistic for rocky planets, but it's far worse for ice and gas giants.

I know that there's been a lot of revision recently on the composition and layers of Jupiter and Saturn, and a lot is based on our ability to model metallic hydrogen about which we know very little. So recent sources would be better than old sources (though I'll take what I can get).

Attempt

Arturo don Juan and A. C. A. C. asked how much this effect matters, so here are some calculations. I'll use Jupiter as my example. Its mass is 1.8986e27 kg and its mean radius is 69911 km, giving it $U_\circ = 2.065 \cdot 10^{36}\text{ J}.$ Using the crude approximation $$ \rho_x=\begin{cases} 25 & \text{ if }x \le .14\\ 6.0533 - 6.3166x & \text{ if }.14<x<.8\\ \frac52(1-x) & \text{ if }x \ge .8 \end{cases} $$ (where 0 is the core and 1 is the edge of the clouds and the density is in g/cm^3) based on the diagram here I compute $$ U \approx 1.775 \cdot 10^{36}\text{ J}. $$

Background

I recently answered a question about the Death Star's destructive ability in Star Wars. My answer relied heavily on the uniform density approximation for scaling, and I'm not sure how much that affected it. I hoped that starting from reasonable baselines would help, but I wasn't really able to do that because I couldn't find any information on the GBE for planets other than Earth. (I could find plenty of exercises asking students to apply the above approximation of course.)

I have to imagine there's something out there, maybe just giving the best-guess density $\rho_x$ at a given depth $x$, from which one can integrate $$ U = \frac{16\pi^2G}{3}\int_0^r x^4\rho_x^2 dx $$ to get the GBE.

Answer 1

The gravitational potential energy of a polytropic sphere of gas (i..e governed by a polytropic equation of state), with a polytropic index $n$ is given by $$ \Omega = -\frac{3}{5-n}\ \frac{GM^2}{R}$$ See for example here.

A constant density sphere would have to have to be incompressible. Since a polytrope has a pressure $p \propto \rho^{1 +1/n}$, this corresponds to $n=0$ and gives you your result for a uniform sphere.

Higher values of $n$ have more centrally condensed profiles and larger binding energies.

The appropriate value of the polytropic index for brown dwarfs and gas giants, where energy transport is convective or where the gas is (non-relativistically) degenerate is $n=3/2$. In which case, your leading coefficient increases from 3/5 to 6/7. i.e. No big deal for an approximate calculation.

Answer 2

This is an answer to the point I brought up in the comments.

To begin, I'd like to point out that the quoted formula for the gravitational binding energy of a spherically symmetric body isn't quite right. You can see in the wiki article that the factor of $3$ in the denominator comes from assuming a uniform mass distribution. Just to have it on record, here's a proper derivation.

To send a layer of mass off the surface of a sphere - its shell - to infinity takes $dU$ energy:

$$dU(r)=\frac{GM_{\text{inside}}(r)\,m_{\text{shell}}(r)}{r}=\frac{GM_{\text{inside}}(r)}{r}(4\pi r^2 \rho(r)\,dr)$$

The mass inside this shell is

$$M_{\text{inside}}(r)=\int_0^r\rho(r') \,(4\pi r'^2\, dr')$$

To get the total binding energy, we send all the shells off to infinity, i.e. we integrate from the core ($r=0$) to the outer radius ($r=R$):

$$U_{\text{binding}}=\int dU = 16\pi^2 G \int_0^R dr_1 \rho(r_1) r_1 \int_0^{r_1} dr_2\, \rho(r_2) r_2^2$$

Maybe with some trick you could reduce that double integral to a single integral, but I can't think of how to do it off the top of my head.

---

If we assume the distribution goes like, say, $\rho(r)=\rho_0 e^{-\alpha r}$, where $\rho_0$ is fixed by the condition that the total mass is $M$, you can carry out the lengthy calculations and find the binding energy as a function of the parameter $\alpha$. Setting $R=M=G=1$, my calculations gave me:

$$U_{\text{binding}}(\alpha)=\frac{\alpha}{4}\frac{\frac{5}{8}-2e^{-\alpha}(\alpha+1)+\frac{1}{8}(2\alpha(\alpha(2\alpha+7)+11)+11)e^{-2\alpha}}{\left(1-(1+\alpha+\frac{1}{2}\alpha^2)e^{-\alpha}\right)^2}$$

Here's a plot of it (courtesy of wolframalpha):

First, as a sanity check, notice that at $\alpha=0$ we recover the uniform-density answer,

$$U_{\text{binding}}(\alpha=0)=\frac{3 M^2 G}{5R}$$

What this plot tells us is that if we have a mass distribution which decays like $e^{-4r}$, the binding energy will change by $\approx (0.8-0.6)/0.6\times 100\%=33\%$. I have specifically just chosen $\alpha=4$ because that gives approximately the mass-density of Jupiter in the link you provided (ignoring the gaseous hydrogen part, whose inclusion would only make the difference more profound).

However, for a rocky planet like earth, the density near the surface is approximately a factor $3/13\approx e^{-1.5}$ less than the density at the core, then $U_{\text{binding}}(\alpha=1.5)=0.649$, which is only $8\%$ different than the uniform-density approximation.

---

Conclusion: The exact mass distribution for a gas giant does matter, because the calculated binding energy of a planet semi-realistic distribution (exponential decay) differs from the uber-naive distribution (uniform) by a little more than $33\%$. However, for a rocky planet whose distribution tapers off much more slowly, the semi-realistic distribution only changes the naive answer about $8\%$.

So, approximating the Earth by a sphere of uniform mass-density isn't so bad, unless you want the answer to better than $8\%$.

---

Next Step: A possible next step would be to do a piecewise distribution, using a different distribution for each "layer". I don't think this would change the answer too much though, but it's a possible direction.