I don't think it is difficult to derive analytically the shape of the Earth. Simply look for the shape of the surfaces of equal potential.
The geometrical symmetry reduces the calculation to a 2-dimensional problem. Assume the rotation axis is vertical. The potential is the sum of the gravitational plus centrifugal:
$\Phi=\Phi_{g}+\Phi_{c}=-\frac{GM_{(x,y)}}{\sqrt{x^2+y^2}}+\frac{\omega^{2}}{2}x^{2}=-\frac{GM_{r}}{r}+\frac{\omega^{2}}{2}r^{2} \cos^{2} l$
The angle $l$ is the same as the latitude, and $M_{r}$ is the mass enclosed by a spherical surface (but please see footnote) at the point, i.e $M_{r}=\frac{4}{3}\pi \rho r^{3}$ by assuming a constant density model. Therefore,
$\Phi= -\frac{4}{3} G \pi \rho r^{2} +\frac{\omega^{2}}{2}r^{2} \cos^{2} l = r^{2}(\frac{\omega^{2}}{2}\cos^{2} l -\frac{4}{3} G \pi \rho)$
Thus, the family of curves of constant (negative) potential $\Phi=-C^{2}$ is:
$ -C^{2} = r^{2}(\frac{\omega^{2}}{2}\cos^{2} l -\frac{4}{3} G \pi \rho) = r^{2}(A^{2} \cos^{2} l -B^{2}) $
Let's go back to rectangular coordinates, to see that this is indeed an ellipse:
$ C^{2} = r^{2}(B^{2} - A^{2} \cos^{2} l) = (x^{2}+y^{2})(B^{2} - A^{2} \frac{x^{2}}{x^{2}+y^{2}}) = (x^{2}+y^{2})B^{2} - A^{2} x^{2}$
$ C^{2} = (B^{2} - A^{2}) x^{2} + B^{2} y^{2} $
For that equation to be an ellipse, $B^{2} - A^{2}$ must be positive. This is natural, otherwise (see how we defined $A$ and $B$) the angular speed $\omega$ would make the centrifugal force stronger than the gravitational force. The semiaxis are then $1/B$ for the vertical direction, and $1/\sqrt{B^{2} - A^{2}}$, i.e. bigger, in the horizontal direction. Note too, that $A=0$ for $\omega = 0$, that is, you recover the spherical shape if there is no rotation.
Thus, an Earth with constant density that rotates as a rigid solid can be approximated by an ellipsoid shape, whose dimension along the rotation axis is smaller.
Additionally, we probably don't need the interior of the Earth to be molten, for the hydrostatic equilibrium assumption to be valid. It could be completely cold and solid and the model still would hold, because at that size scales, relative small deviations of matter distribution from the constant potential surfaces give rise to enormous shear stress that rocks, no matter how hard and solid, cannot resist. That is why the liquid model is a valid approximation (but I have not done any numbers on this).
NOTE: We have assumed that any point belongs to a spherical surface that is completely full of matter, therefore the potential gravitational energy is the same as if all matter inside that sphere were located at the Earth centre. If the Earth were much more flattened, this approximation would not be valid.
g
that is measured at sea level on poles or equator. $\endgroup$