A simple proof that under Newtonian gravity rotating massive bodies are ellipsoids?
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)}}. 521
$$ 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? 521