Ok, here is the promised follow-up post on the case of a spherically symmetric mass of uniform density. I should say at the outset that digging into this has made me realize that the interpretation of ##M_0 - M## as "gravitational binding energy" is at best heuristic; the reasoning involved has a lot of caveats in it, and I have updated the Insights post accordingly (I have also updated it with some other corrections arising from what follows). However, this also gives me some good material for the next follow-up post in the series, so it all works out in the end.
We are using units in which ##G = c = 1##. As I said, this case has a known exact solution for the metric, which is as follows:
$$
ds^2 = - \left[ \frac{3}{2} \sqrt{1 - \frac{8}{3} \pi \rho R^2} - \frac{1}{2} \sqrt{1 - \frac{8}{3} \pi \rho r^2} \right]^2 dt^2 + \frac{1}{1 - \frac{8}{3} \pi \rho r^2} dr^2 + r^2 d\Omega^2
$$
where ##\rho## is the density (and is constant) and ##R## is the radial coordinate of the object's surface (i.e., for ##r > R## there is vacuum). For convenience, we define a constant ##k = \frac{8}{3} \pi \rho##; the metric then becomes
$$
ds^2 = - \left[ \frac{3}{2} \sqrt{1 - k R^2} - \frac{1}{2} \sqrt{1 - k r^2 } \right]^2 dt^2 + \frac{1}{1 - k r^2} dr^2 + r^2 d\Omega^2
$$
The Killing vector field ##\xi^a## is just ##\partial / \partial t## in these coordinates, and its norm is ##\sqrt{\xi^a \xi_a} = \sqrt{g_{ab} \xi^a \xi^b} = \sqrt{g_{tt}}##, so
$$
\sqrt{\xi^a \xi_a} = \left[ \frac{3}{2} \sqrt{1 - k R^2} - \frac{1}{2} \sqrt{1 - k r^2 } \right]
$$
The 4-velocity ##u^a## is equal to ##\xi^a / \sqrt{\xi^a \xi_a}##, so its only nonzero component is
$$
u^0 = \frac{1}{\left[ \frac{3}{2} \sqrt{1 - k R^2} - \frac{1}{2} \sqrt{1 - k r^2 } \right]}
$$
However, we won't actually need this in this particular case, as we'll see in a moment.
Next, we have the stress-energy tensor ##T_{ab}## and its trace ##T##; the quantity that appears in the mass integrals is ##\left( 2 T_{ab} - g_{ab}T \right) u^a u^b = \left( \rho + 3 p \right)##, where ##p## is the pressure. (Note that this is really a way of stating the physical definition of ##\rho## and ##p##, because we are contracting the stress-energy tensor with the 4-velocity--this is why we don't need to actually know the expression for the 4-velocity components.) However, even though the density ##\rho## is constant, the pressure ##p## is not; it increases as ##r## decreases, from ##p = 0## at the surface to some positive value ##p_c## at the center ##r = 0##. So we need the equation for ##p(r)##, which turns out to be
$$
p = \rho \frac{\sqrt{1 - k r^2} - \sqrt{1 - k R^2}}{3 \sqrt{1 - k R^2} - \sqrt{1 - k r^2}}
$$
Finally, we need to look at the volume element ##dV##. In flat spacetime, this volume element, after integrating over the angular directions, would simply give ##4 \pi r^2 dr##, i.e., the area of a 2-sphere at radius ##r## times the radial differential. However, in curved spacetime, this coordinate volume element is not actually the physical volume element; there is an extra factor of ##\sqrt{g_{rr}}##, to account for curvature. (I should have made this clearer in previous posts, since it has implications for the physical meaning of ##M_0##; see further comments below.) So we have
$$
\int dV = \int 4 \pi r^2 \frac{1}{\sqrt{1 - k r^2}} dr
$$
We now have all we need to evaluate the integrals. We will do the integral for ##M## first, for reasons which will shortly become apparent; it is
$$
M = \int \sqrt{\xi^a \xi_a} \left( 2 T_{ab} - g_{ab} T \right) u^a u^b dV = \int_0^R \left[ \frac{3}{2} \sqrt{1 - k R^2} - \frac{1}{2} \sqrt{1 - k r^2 } \right] \left( \rho + 3 p \right) 4 \pi r^2 \frac{1}{\sqrt{1 - k r^2}} dr
$$
Substituting for ##p## gives
$$
M = \int_0^R \left[ \frac{3}{2} \sqrt{1 - k R^2} - \frac{1}{2} \sqrt{1 - k r^2 } \right] \rho \left[ 1 + 3 \frac{\sqrt{1 - k r^2} - \sqrt{1 - k R^2}}{3 \sqrt{1 - k R^2} - \sqrt{1 - k r^2}} \right] 4 \pi r^2 \frac{1}{\sqrt{1 - k r^2}} dr
$$
The messy expressions involving the square roots turn out to all cancel, leaving
$$
M = \int_0^R 4 \pi \rho r^2 dr = \frac{4}{3} \pi \rho R^3
$$
Notice that this is what we would have expected in Newtonian physics! That means I actually told a bit of fib in an earlier post, when I said ##M## was "different" in GR vs. NP; I had forgotten that in GR, pressure gravitates. We'll discuss that further below. However, it was only a bit of a fib, because even though this result formally looks the same as the Newtonian result, the way we arrived at it makes clear that we did not just "naively" integrate ##\rho## over the Euclidean volume of the object. We integrated a much more complicated expression, involving effects of both pressure (in addition to density) and spacetime curvature, in which those effects just happened to cancel in the right way. Again, more on this below.
The fact that ##M## came out this way does illustrate, though, what I said in an earlier post about the difference between ##M_0## and ##M## being negligible if we really adopt the weak field, slow motion limit. In that limit, if we do indeed ignore spacetime curvature, and if we also ignore pressure, because it will be negligible compared to energy density in this limit, then the integral for ##M_0## will in fact look exactly like the above--and in fact it will also be the same as the "naive" integral that we started out with in the Insights post! This is because ##T_{ab} u^a u^b = \rho##--this is the physical definition of ##\rho##--and so ##M_{naive}## from the post also looks exactly the same as the integral for ##M## above. This is one of the key caveats in trying to interpret what ##M_0## means, physically.
If we don't ignore pressure, but we do leave out the effects of spacetime curvature, we can do the integral for ##M_0## and expect to get something different (and hopefully larger) than what we got for ##M## above; but we have to consider, first, what volume element ##dV## we should use for this case. If we are assuming flat spacetime (or at least assuming that the effects of space curvature are negligible because the field is weak enough), then we should not have the extra factor of ##\sqrt{g_{rr}}## in the volume element. (Again, we'll discuss this further below.) So we have:
$$
M_0 = \int \left( 2 T_{ab} - g_{ab} T \right) u^a u^b dV = \int_0^R \left( \rho + 3 p \right) 4 \pi r^2 dr = \int_0^R \rho \left[ 1 + 3 \frac{\sqrt{1 - k r^2} - \sqrt{1 - k R^2}}{3 \sqrt{1 - k R^2} - \sqrt{1 - k r^2}} \right] 4 \pi r^2 dr
$$
Here we don't have quite the same nice cancellation, but we can still simplify quite a bit:
$$
M_0 = \int_0^R \rho \left[ \frac{2 \sqrt{1 - k r^2}}{3 \sqrt{1 - k R^2}} \frac{1}{1 - \frac{\sqrt{1 - k r^2}}{3 \sqrt{1 - k R^2}}} \right] 4 \pi r^2 dr
$$
This is still quite messy, but we can use the fact that, for weak fields and slow motion, ##k R^2 < k r^2 << 1##, so we can expand out the denominators and square roots and drop terms of second and higher order in ##k## to obtain, after some algebra:
$$
M_0 = M + \frac{4}{27} \frac{M^2}{R}
$$
I'm not entirely sure all of the algebra is correct here, but we expect the leading term to be ##M##, and the correction term is based on how the various square roots and denominators expand out and then get integrated over ##r##. This correction, of course, is quite a bit smaller than the expected Newtonian value of ##3 M^2 / 5 R##, which at least shows, as I said in a previous post, that we should not expect the Newtonian values to always be good estimates of what a GR calculation will show, even in the weak field, slow motion limit.
However, let's now discuss what, if anything, the ##M_0## calculation we just did actually means, physically. On its face, the calculation is certainly "wrong", in the sense that we ignored spacetime curvature but did not ignore pressure, even though the effects of both are expected to be of the same order (we know that because those effects canceled exactly in the integral for ##M## above). But, as we saw, if we ignore both spacetime curvature and pressure, we get the same answer as ##M## above--i.e., we don't even see any "gravitational binding energy" at this order of approximation. And if we include both spacetime curvature and pressure, we get ##M## by the route we got it above, because their effects cancel.
But let's now consider the alternative interpretation of what ##M_0## might mean, as I described in a previous post. Imagine that we have a bound object, with mass ##M## as calculated by the correct integral above, and we now "disassemble" it into its constituent particles, moving each particle out to infinity (meaning, in practice, far enough away from all other particles that gravity between them is negligible). We will obviously have to add energy to the system to do this; and at the end of it, we will have a system composed of a cloud of widely separated particles, each of which has negligible self-gravity. This cloud will have some density ##\bar{\rho}##, and will be contained within a spherically symmetric region of some radius ##\bar{R}##, and its total mass should be given by
$$
\bar{M} = \int_0^\bar{R} 4 \pi \bar{\rho} r^2 dr = \frac{4}{3} \pi \bar{\rho} \bar{R}^3
$$
This should be true because the cloud is so dispersed that its spacetime curvature should really be negligible, and its pressure should also really be negligible since the particles are too far apart to interact with each other, so we really can use the "naive" integral in this case to get the right answer.
Then, if ##E## is the total energy we had to add to the original object to "disassemble" it into the widely dispersed cloud, we should have
$$
E = \bar{M} - M = \frac{4}{3} \pi \left( \bar{\rho} \bar{R}^3 - \rho R^3 \right)
$$
In other words, on this intrepretation, the difference between ##M_0## (which we are calling ##\bar{M}## now) and ##M## arises from the fact that, when we "disassemble" a bound object, the product of its density and its radius cubed increases--the density doesn't decrease as much as the cube of the radius increases. But even considering this interpretation requires relativistic energy-mass equivalence; as I've said before, in Newtonian physics, the energy ##E## being added does not change the total mass, and a Newtonian calculation would predict ##\bar{M} = M##, i.e., it would predict that, as an object is "disassembled", its density should decrease in exact proportion to the increase in the cube of its radius. In GR, however, this doesn't happen; adding the energy ##E## really does add mass to the system, and that added mass shows up in the slower decrease of the density as the cube of the radius increases, so ##\bar{M} > M##.
This has been a long post, but hopefully it has helped to clarify what is going on in this scenario.
@exponent137, your original question has sparked a very good discussion, and that's the mark of a good question!
