Shapes and gravitational fields of rotating two-layer Maclaurin   ellipsoids: Application to planets and satellites
  The exact solution for the shape and gravitational field of a rotating two-layer Maclaurin ellipsoid of revolution is compared with predictions of the theory of figures up to third order in the small rotational parameter of the theory of figures. An explicit formula is derived for the external gravitational coefficient $J_2$ of the exact solution. A new approach to the evaluation of the theory of figures based on numerical integration of ordinary differential equations is presented. The classical Radau-Darwin formula is found not to be valid for the rotational parameter \epsilon_2 = \Omega^2/(2\pi G\rho_2) >= 0.17 since the formula then predicts a surface eccentricity that is smaller than the eccentricity of the core-envelope boundary. Interface eccentricity must be smaller than surface eccentricity. In the formula for $\epsilon_2$, $\Omega$ is the angular velocity of the two-layer body, $\rho_2$ is the density of the outer layer, and G is the gravitational constant. For an envelope density of 3000 kg m^(-3) the failure of the Radau-Darwin formula corresponds to a rotation period of about 3 hr. Application of the exact solution and the theory of figures is made to models of Earth, Mars, Uranus, and Neptune. The two-layer model with constant densities in the layers can provide realistic approximations to terrestrial planets and icy outer planet satellites. The two-layer model needs to be generalized to allow for a continuous envelope (outer layer) radial density profile in order to realistically model a gas or ice giant planet. 
Gerald Schubert, John D. Anderson, Keke Zhang, Dali Kong, and Ravit   Helled
astro-ph.EP