How to find the density of a planet and its core taking into account the gravitational compression in them?

Question

Consider a planet $P$ of radius $r$ larger than Earth, whose mass is $r^3$. By the square–cube law, this planet will have the same density as the Earth... assuming that the gravitational compression on it is the same as on our planet, but this is wrong, because the planet will have a greater gravity, therefore, the gravitational compression in it will be more intense, making it smaller and denser than expected (is this reasoning correct?). But what matters to me is not only the planet's density, but also that of its core, which would also increase.

The Earth's core has a density of 11,000 kg/m³, greater than would be expected if it were calculated considering it only as a mixture of two substances (Fe and Ni, to simplify), without taking into account the gravitational compression. How can I calculate the density of the core of $P$ (and of itself) taking into account the gravitational compression? Is it possible to do it without much trouble?

Answer 1

The technique to do this is similar to that used in constructing stellar models. You know some of the properties of your object - in this case, it look like the mass and radius. You want to figure out the internal structure of the planet, including the central density and pressure, as well as the density and pressure profiles as functions of radius. The best way to do this is as follows:

1. Determine a relationship between density and pressure - what we call an equation of state. We write this as $\rho = f(P)$.
2. Make a guess as to the central pressure - and therefore the central density, using the above relation.
3. Numerically integrate the differential equations determining the structure of the planet. Start in the core, and go outwards.
4. Once your equations indicate that you've reached the surface of the planet, see what the predicted mass and radius are.
5. Given the result of step 4, adjust your guesses for the central pressure and density, and repeat steps 3 and 4.

The two equations you need, besides the equation of state, are the equation of hydrostatic equilibrium $(1)$ and the equation of mass continuity $(2)$: $$\frac{dP}{dr} = -\frac{GM_r\rho}{r^2}\tag{1}$$ $$\frac{dM}{dr} = 4\pi r^2\rho\tag{2}$$ Here, $P$ is pressure, $\rho$ is density, $r$ is radius, $G$ is the gravitational constant, and $M_r$ is the mass contained within a radius $r$: $$M_r = \int_0^r 4\pi r'^2\rho(r')dr'$$ Let's say that at a given radius $r_i$, you know the pressure, density, and mass encapsulated within that radius - $P_i$, $\rho_i$, and $M_{r,i}$. Assume that your integration uses some step size $\Delta r$. Then the values at the radius $r_{i+1}$ are $$r_{i+1} = r_i + \Delta r$$ $$P_{i+1} = P_i + \Delta P = p_i + \frac{dP}{dr}\Delta r = P_i - \frac{GM_{r,i}\rho_i}{r_i^2}\Delta r$$ $$\rho_{i+1} = f(P_{i+1})$$ $$M_{r,i+1} = M_i + \Delta M = M_i + \frac{dM}{dr}\Delta r = M_i + 4\pi r_i^2\rho_i\Delta r$$ Continue doing this until you reach the surface of the planet, i.e. at $P = 0$, then repeat steps 3 to 4 with a new guess until you're satisfied with the model's predicted mass and radius.

I implemented this procedure in my answer to another question, although I only had to do steps 3 to 4 once - I had a reasonable central pressure all figured out. I used an equation of state from Seager et al. 2008, designed for rocky planets.