A simple proof that under Newtonian gravity rotating massive bodies are ellipsoids?

Question

Here is my attempt at deriving the shape of an idealized rotating massive body under Newtonian gravity, assuming that the gravity force points towards the center of mass and shape of the body is stabilized. Physically one might imagine an object with very heavy core localized in very small volume at the center, surrounded by very light liquid.

A point on the surface moves according to $\ell=\omega t$, where $t$ is time, $\ell$ is longitudal angle and $\omega$ the constant angular velocity. Thus in cartesian coordinates $$ {\mathbf r}(h,t)=R(h)(\cos(\omega t)\cos(h),\sin(\omega t)\cos(h),\sin(h)) $$ where $h$ is the latitude (angle) and $R(h)$ is distance to the center of mass (which by assumption does not depend on time, hence also on the longitude since the latter depends linearly on time).

Net force on a mass $m$ at this point is thus $$ {\mathbf F}=m{\mathbf r}_{tt}=-\omega^2mR(h)\cos(h)(\cos(\omega t),\sin(\omega t),0) $$

If that mass does not participate in any other motion except the overall rotation, then this net force must be equal to the sum of the gravity force and the normal force at the (frictionless) surface of the planet.

The gravity force (with gravitational constant $G$ and planet mass $M$) is $$ {\mathbf F}_g=-\frac{GMm}{R(h)^2}(\cos(\omega t)\cos(h),\sin(\omega t)\cos(h),\sin(h)) $$ and we have to write down the condition that $\mathbf F-{\mathbf F}_g$ points normally to the surface.

One normal vector is \begin{multline*} {\mathbf r}_\ell\times{\mathbf r}_h=\\R\cos(h)((R\cos(h)+R'\sin(h))\cos(\omega t),(R\cos(h)+R'\sin(h))\sin(\omega t),R\sin(h)-R'\cos(h)) \end{multline*} (I shortened $R(h)$ to $R$ and $\frac{dR(h)}{dh}$ to $R'$ just to fit the expression on one line.)

The condition that $\mathbf F-{\mathbf F}_g$ is parallel to it gives $$ \frac{R\cos(h)+R'\sin(h)}{\frac{GM}{R^2}-\omega^2R}=\frac{R\sin(h)-R'\cos(h)}{\frac{GM}{R^2}\tan(h)}; $$ solving this for $R'$ we arrive at the differential equation $$ R'=\frac1{\frac1{R\tan(h)}-\frac{2GM}{\omega^2}\frac1{R^4\sin(2h)}}. $$ Solutions of the latter are determined by $$ (R\cos(h))^2=\frac{2GM}{\omega^2}(C-\frac1R) $$ with arbitrary constant $C$. Switching to $x=R\cos(h)$, $y=R\sin(h)$ and denoting $\frac{2GM}{\omega^2}=k$, $kC=a$ we obtain $$ y=\pm\sqrt{\left(\frac k{a-x^2}\right)^2-x^2}. $$

And now I am stuck since this is most certainly not an ellipsoid.

Where is my error? Can all this mess be simplified?

Answer 1

As you mention, you are seeking to set up calculation for the idealized case.

As we know: actual celestial bodies do not have a homogeneous density; you are aware that treating the rotating celestial body as having a homogenous density is unphysical, but without that simplification the problem would be untractable.

But you also state another assumption, one that is actually invalid. It is only in the case of a perfectly spherical celestial body that all over the surface the true gravity points towards the geometric center of the celestial body.

I use the expressions 'true gravity' and 'effective gravity' for the following distinction: effective gravity is the resultant of true gravity and the requirement to provide centripetal force to make objects co-rotate with the celestial body.

In the case of the Earth there are two contributions to the fact that on the Equator the effective gravitational acceleration is smaller than at the poles.

• The acceleration that is required to circumnavigate the Earth's axis along the equator at one revolution per sidereal day is 0.0339 $\tfrac{m}{s^2}$. That goes at the expense of true gravity.

• At the Equator you are about 21 kilometers further away from the geometric center of the Earth than at the poles. The corresponding difference in gravitational acceleration is 0.0178 $\tfrac{m}{s^2}$

So: the distribution of those two contributions is about 70/30

When you are at the Equator the center of gravitational attraction (true gravity) on you does not coincide with the geometric center of the Earth.

It's off by less than 10 kilometers, but its very much not negligable because it is on the same order of size as the deviation that you want to compute: the oblateness of the celestial body.

This extends to all latitudes; tor every latitude the center of gravitational attraction is located a bit differently.

The starting point to solve for the shape of the rotating celestial body has to be that the rotating body is in hydrostatic equilibrium. As you mention: everywhere along the surface the local effective gravity is perpendicular to the local horizontal.

You can set up an expression that gives for every point along the surface the required centripetal acceleration.

For a rotating celestial body: the equatorial radius is larger than the polar radius. In effect the surface from the Equator to the pole is a downhill slope.

I expect it is possible to set up an expression for the magnitude of that slope as a function of latitude. Solving for that slope as a function of latitude will, I expect, give the shape of the rotating celestial body.

Answer 2

OK so I now managed to slightly simplify it and I confirm that it is not an ellipsoid.

Presuming axial symmetry it suffices to consider only two coordinates: distance to the axis $x$ and distance to the equatorial plane $z$.

Net force on a point of mass $m$ at rest on the surface with coordinates $(x,z)$ is the centripetal force $$ -\omega^2m(x,0) $$ which must be the sum of a normal force and the gravity force $$ -\frac{GMm}{R^2}(x/R,z/R) $$ where $R^2=x^2+z^2$. Thus $$ -\omega^2m(x,0)-(-\frac{GMm}{R^2}(x/R,z/R)) $$ must be perpendicular to the surface. Tangent to the surface is $(dx,dz)$, so we obtain $$ (-\omega^2mx+\frac{GMm}{R^2}x/R)dx+(\frac{GMm}{R^2}z/R)dz=0, $$ i. e. $$ \frac{GM}{R^3}(xdx+zdz)=\omega^2xdx. $$ Using $xdx+zdz=RdR$ this is equivalent to $$ (dR)/R^2=\varepsilon d(x^2) $$ where $$ \varepsilon=\frac{\omega^2}{2GM} (\approx 5\times10^{-29} \text{ for Earth}). $$ As $(dR)/R^2=-d(1/R)$, we obtain $$ c-1/R=\varepsilon x^2 $$ for some constant $c$. When $\varepsilon=0$ we must have $c=1/R_0$ where $R_0$ is the radius in the absence of rotation, when the surface is spherical.

We finally obtain $$ x^2+z^2=\frac{R_0^2}{(1-\varepsilon x^2)^2}. $$ Here is a plot for several values of $\varepsilon$ small compared to $R_0$, which shows that this is definitely not an ellipse.