Why is the Earth so fat?

Question

I made a naive calculation of the height of Earth's equatorial bulge and found that it should be about 10km. The true height is about 20km. My question is: why is there this discrepancy?

The calculation I did was to imagine placing a ball down on the spinning Earth. Wherever I place it, it shouldn't move.

The gravitational potential per unit mass of the ball is $hg$, with $h$ the height above the pole-to-center distance of the Earth (call that $R$) and $g$ gravitational acceleration.

Gravity wants to pull the ball towards the poles, away from the bulge. It is balanced by the centrifugal force, which has a potential $-\omega^2R^2\sin^2\theta/2$ per unit mass, with $\omega$ Earth's angular velocity and $\theta$ the angle from the north pole. This comes taking what in an inertial frame would be the ball's kinetic energy and making it a potential in the accelerating frame.

If the ball doesn't move, this potential must be constant, so

$$U = hg - \frac{(\omega R \sin\theta)^2}{2} = \textrm{const}.$$

We might as well let the constant be zero and write

$$h = \frac{(\omega R \sin\theta)^2}{2g}.$$

For the Earth,

$$R = 6.4*10^6m$$

$$\omega = \frac{2\pi}{24\ \textrm{hours}}$$

$$g = 9.8\ m/s^2.$$

This gives 10.8 km when $\theta = \pi/2$, so the equatorial bulge should be roughly that large.

According to Wikipedia, the Earth is 42.72 km wider in diameter at the equator than pole-to-pole, meaning the bulge is about twice as large as I expected. (Wikipedia cites diameter; I estimated radius.)

Where is the extra bulge coming from? My simple calculation uses $g$ and $R$ as constants, but neither varies more than a percent or so. It's true the Earth does not have uniform density, but it's not clear to me how this should affect the calculation, so long as the density distribution is still spherically-symmetric (or nearly so).

(Wikipedia also includes an expression, without derivation, that agrees with mine.)

Answer 1

The error is that you assume that the density distribution is "nearly spherically symmetric". It's far enough from spherical symmetry if you want to calculate first-order subleading effects such as the equatorial bulge. If your goal is to compute the deviations of the sea level away from the spherical symmetry (to the first order), it is inconsistent to neglect equally large, first-order corrections to the spherical symmetry on the other side - the source of gravity. In other words, the term $hg$ in your potential is wrong.

Just imagine that the Earth is an ellipsoid with an equatorial bulge, it's not spinning, and there's no water on the surface. What would be the potential on the surface or the potential at a fixed distance from the center of the ellipsoid? You have de facto assumed that in this case, it would be $-GM/R+h(\theta)g$ where $R$ is the fixed Earth's radius (of a spherical matter distribution) and $R+h(\theta)$ is the actual distance of the probe from the origin (center of Earth). However, by this Ansatz, you have only acknowledged the variable distance of the probe from a spherically symmetric source of gravity: you have still neglected the bulge's contribution to the non-sphericity of the gravitational field.

If you include the non-spherically-symmetric correction to the gravitational field of the Earth, $hg$ will approximately change to $hg-hg/2=hg/2$, and correspondingly, the required bulge $\Delta h$ will have to be doubled to compensate for the rotational potential. A heuristic explanation of the factor of $1/2$ is that the true potential above an ellipsoid depends on "something in between" the distance from the center of mass and the distance from the surface. In other words, a "constant potential surface" around an ellipsoidal source of matter is "exactly in between" the actual surface of the ellipsoid and the spherical $R={\rm const}$ surface.

I will try to add more accurate formulae for the gravitational field of the ellipsoid in an updated version of this answer.

Update: gravitational field of an ellipsoid

I have numerically verified that the gravitational field of the ellipsoid has exactly the halving effect I sketched above, using a Monte Carlo Mathematica code - to avoid double integrals which might be calculable analytically but I just found it annoying so far.

I took millions of random points inside a prolate ellipsoid with "radii" $(r_x,r_y,r_z)=(0.9,0.9,1.0)$; note that the difference between the two radii is $0.1$. The average value of $1/r$, the inverse distance between the random point of the ellipsoid and a chosen point above the ellipsoid, is $0.05=0.1/2$ smaller if the chosen point is above the equator than if it is above a pole, assuming that the distance from the origin is the same for both chosen points.

Code:

{xt, yt, zt} = {1.1, 0, 0};

runs = 200000;
totalRinverse = 0;
total = 0;

For[i = 1, i < runs, i++,
 x = RandomReal[]*2 - 1;
 y = RandomReal[]*2 - 1;
 z = RandomReal[]*2 - 1;
 inside = x^2/0.81 + y^2/0.81 + z^2 < 1;
 total = If[inside, total + 1, total];
 totalRinverse = 
  totalRinverse + 
   If[inside, 1/Sqrt[(x - xt)^2 + (y - yt)^2 + (z - zt)^2], 0];
]

res1 = N[total/runs / (4 Pi/3/8)]
res2 = N[totalRinverse/runs / (4 Pi/3/8)]
res2/res1

Description

Use the Mathematica code above: its goal is to calculate a single purely numerical constant because the proportionality of the non-sphericity of the gravitational field to the bulge; mass; Newton's constant is self-evident. The final number that is printed by the code is the average value of $1/r$. If {1.1, 0, 0} is chosen instead of {0, 0, 1.1} at the beginning, the program generates 0.89 instead of 0.94. That proves that the gravitational potential of the ellipsoid behaves as $-GM/R - hg/2$ at distance $R$ from the origin where $h$ is the local height of the surface relatively to the idealized spherical surface.

In the code above, I chose the ellipsoid with radii (0.9, 0.9, 1) which is a prolate spheroid (long, stick-like), unlike the Earth which is close to an oblate spheroid (flat, disk-like). So don't be confused by some signs - they work out OK.

Bonus from Isaac

Mariano C. has pointed out the following solution by a rather well-known author:

http://books.google.com/books?id=ySYULc7VEwsC&lpg=PP1&dq=principia%20mathematica&pg=PA424#v=onepage&q&f=false

Answer 2

I) Flatness. Here, we would like to calculate analytically Lubos Motl's solution to the first order in the flatness parameter $f$,

$$0~<f~:=~1-\frac{b}{a}~\approx~\frac{a}{b}-1 ~\ll~ 1, \tag{1}$$

where $a$ and $b$ are the equatorial and polar radius of the Earth, respectively, and $a>b$. (The $\approx$ symbol will, from now on, mean equality up to higher-order terms in $f$.) We assume that Earth is a massive oblate ellipsoid

$$ \left(\frac{x}{a}\right)^2 + \left(\frac{y}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~\leq~ 1 \tag{2} $$

with uniform density $\rho$ and volume

$$V~=~\frac{4}{3}\pi a^2 b~\approx~\frac{4}{3}\pi b^3(1+2f).\tag{3}$$

The eccentricity is

$$e~=~\sqrt{(2-f)f}~\approx~ \sqrt{2f}~\ll ~1.\tag{4}$$

II) Quadrupole. We will assume that gravity is Newtonian. What we need to calculate is the quadrupole moment contribution to the gravitational potential

$$ U_4 ~= ~\frac{x^i Q_{ij}x^j}{r^5},\tag{5} $$

For symmetry reasons, one of the principal directions of the quadrupole moment $Q_{ij}=Q_{ji}$ must be along the polar $z$-axis, and the two other principal directions must have the same eigenvalues and be in the equatorial $xy$ plane. Thus, $U_4$ must be of the form

$$ U_4 ~=~ \frac{Q_a s^2+Q_b c^2}{r^3}, \tag{6}$$

where $Q_a$ and $Q_b$ are the equatorial and polar eigenvalues, respectively; where we have introduced the short-hand notation $s\equiv\sin(\theta)$ and $c\equiv\cos(\theta)$; and where $\theta\in[0,\pi]$ is the polar angle. Since the quadrupole moment $Q_{ij}$ cannot contribute to Gauss' law, we must demand that the Laplacian $\nabla^2 U_4 = 0$ vanishes, which leads to that the polar eigenvalue is minus two times the equatorial eigenvalue, $Q_b=-2Q_a$. In other words,

$$ U_4 ~\sim~ \frac{3c^2-1}{r^3}~=~\frac{2-3s^2}{r^3}.\tag{7} $$

It is, therefore, enough to calculate the gravitational potential at a point $P=(0,0,r)$ on the polar $z$-axis, where we have azimuthal $U(1)$ symmetry $\varphi\to\varphi+\varphi_{0}$.

III) Ring. Let us calculate the contribution to the gravitational potential energy from a ring parallel to the $xy$ plane and with center $(0,0,z)$ on the $z$-axis. Because of the $U(1)$ azimuthal symmetry, we may focus on a point on the ring with $y=0$ and $x\geq 0$, and which lie in the $xz$ plane. Let this point have 2D polar coordinates $(R,\theta)$. In other words, the point has 3D coordinates $(x,y,z)=(Rs,0,Rc)$. From the equation of an ellipse

$$ \left(\frac{x}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~=~ 1, \tag{8} $$

we obtain, after some elementary algebra,

$$ R ~\approx~ b (1 + fs^2) \quad \left\{\begin{array}{rcl}~=~ b &\mathrm{for}& \theta~=~0, \\ ~\approx~ a &\mathrm{for}& \theta~=~\frac{\pi}{2} . \end{array}\right.\tag{9}$$

Next, the distance $\ell$ from the ring to the point $P$ is given by the cosine relation

$$ \ell^2 ~=~ r^2 +R^2 - 2rRc ~\approx~ r^2 +b^2(1 + 2fs^2) - 2rb(1 + fs^2)c$$ $$~=~(r^2 +b^2 -2rbc) +2fbs^2(b-rc).\tag{10}$$

The "outer" surface area of the ring is

$$ \mathrm d A ~=~ 2\pi Rs \cdot R \mathrm d \theta ~\approx~ 2\pi b^2(1 + 2fs^2) (-\mathrm d c), \tag{11}$$

with infinitesimal "radial" thickness

$$ \mathrm dR ~\approx~ (1 + fs^2) \mathrm db.\tag{12}$$

Thus, the volume of the ring is

$$\mathrm d^2V ~=~ \mathrm dA \ \mathrm dR ~\approx~ 2\pi b^2 \mathrm db(1 + 3fs^2) (-\mathrm dc). \tag{13}$$

(Eventually, we want to integrate over the polar angle $\theta$ from $0$ to $\pi$. This corresponds to integration over $c\equiv\cos(\theta)$ from $1$ to $-1$ in the negative direction. Therefore, $\mathrm dc<0$ is negative.)

IV) Potential. The contribution from the ring to the gravitational potential energy $\mathrm d^2 U$ at the point $P$ is

$$ d^2 U ~=~ -G \frac{\mathrm d^2 M}{\ell} ~=~ - G\rho \frac{\mathrm d^2 V}{\ell} ~\approx~ G\rho\frac{2\pi b^2 \mathrm db(1 + 3fs^2) \mathrm d c}{\sqrt{(r^2 +b^2 -2rbc) +2fbs^2(b-rc)}} $$ $$~\approx~ G\rho2\pi b^2 \mathrm d b\ \mathrm d c\left(\frac{1 + 3f(1-c^2)}{\sqrt{r^2 +b^2 -2rbc}}- \frac{1}{2} \frac{2fb(1-c^2)(b-rc)}{(r^2 +b^2 -2rbc)^{\frac{3}{2}}}\right).\tag{14}$$

Integration over cosine $c$ of the polar angle $\theta\in[0,\pi]$ to form a thin shell yields (with the help of the MAPLE program)

$$ \mathrm d U ~=~ \int_{c=1}^{c=-1} \mathrm d^2 U ~\approx~ -G\rho\frac{4\pi b^2 \mathrm d b}{r}\left( (1+3f) -f \left( 1+ \frac{2b^2}{5r^2}\right) - \frac{4fb^2}{15r^2}\right) $$ $$ ~=~ -G\rho\frac{4\pi b^2 \mathrm d b}{r}\left( (1+2f) - \frac{2fb^2}{3r^2}\right).\tag{15} $$

Next, we integrate over $b$ to get the polar potential of a massive ellipsoid

$$ U ~=~ \int^{b} \mathrm d U ~\approx~ -G\rho\frac{4\pi b^3}{3r}\left( (1+2f) - \frac{2fb^2}{5r^2}\right) ~=~U_1+U_{4b}.\tag{16}$$

V) Monopole & quadrupole. The first term $U_1$ is, of course, precisely the monopole potential

$$U_1(r)~=~-\frac{GM}{r},\tag{17}$$

and the second term $U_{4b}$ is the polar quadrupole potential

$$U_{4b}~\approx~-\frac{2fb^2}{5r^2}U_1.\tag{18}$$

Thus, we know that the gravitational quadrupole potential in an arbitrary point (not necessarily on the $z$-axis) is

$$U_{4}(r)~\approx~-\frac{(2-3s^2)fR^2}{5r^2}U_1(r),\tag{19}$$

and the full gravitational potential is

$$U(r) ~\approx~ U_1(r)+U_4(r) ~\approx~ -\frac{GM}{r}\left(1 - \frac{(2-3s^2)fR^2}{5r^2}\right). \tag{20}$$

VI) Surface. From now on, let us only consider points $P$ with

$$r~=~R~\approx~ b (1 + fs^2)\tag{21}$$

on the surface of the ellipsoid-shaped Earth. Then, the monopole potential in $P$ becomes

$$ U_1(R)~ \approx~ -(1- fs^2)\frac{GM}{b}, \tag{22}$$

while the quadrupole potential in $P$ becomes

$$ U_4(R) ~\approx~ \frac{f(2-3s^2)}{5}\frac{GM}{b} \tag{23}$$

so that the full gravitational potential in $P$ becomes

$$ U(R) ~\approx~ U_1(R) + U_4(R) ~\approx~ - \left(1-\frac{2f(1+s^2)}{5}\right)\frac{GM}{b}.\tag{24}$$

VII) Discussion. Let us use the North Pole as a reference point, i.e., subtract the gravitational potential $U_b(b)$ at the North Pole. Then, the gravitational potential in $P$ becomes

$$ U(R)- U_b(b) ~\approx~ \frac{2fs^2}{5} \frac{GM}{b}~\approx~ \frac{2}{5}g hs^2,\tag{25} $$

where

$$h~:~=~a-b~=~fa~ \approx~ fb ~>~0\tag{26}$$

is the difference in the equatorial and polar radii and

$$g~=~ |\vec{\nabla} U(R)| ~=~ \left. \frac{\partial U(r)}{\partial r}\right|_{r=R}+{\cal O}(f) ~=~ \frac{GM}{R^2}+{\cal O}(f).\tag{27}$$

We should now add the centrifugal potential

$$ V ~=~ -\frac{(\omega Rs)^2}{2}. \tag{28}$$

To the order in $f\ll 1$ that we are working, we see that the total potential $U(R)+V$ is constant (independent of the surface point $P$) if

$$\frac{2}{5}gh~=~\frac{(\omega R)^2}{2}. \tag{29}$$

Conclusion: We find a factor $\color{Red}{\dfrac{2}{5}}$ in difference from Mark Eichenlaub's original monopole argument.

---

Update: Half a year after this answer was posted on Phys.SE, on November 18th, 2011‎, the Wikipedia page changed its listed mathematical expression for the flatness parameter into $\dfrac{h}{R}\approx f\approx\dfrac{5}{4}\dfrac{\omega^2R}{g}$, and is now in full agreement with this answer.

Answer 3

There was some doubt about Lubos' answer (which I've accepted), so this is just a verification.

I copied the method Lubos described and found the potential difference for an ellipsoid with different eccentricities. Sure enough, for an oblate spheroid, if you make the center-equator distance a fraction $e$ larger than the center-pole distance, the potential is roughly a fraction $e/2$ smaller at the equator.

To do the whole problem, we'd have to take into account the varying density of the Earth, but as a rough estimate this seems to do the job.

For example, for the oblate spheroid

$$\left(\frac{x}{1.01}\right)^2 + \left(\frac{y}{1.01}\right)^2 + z^2 < 1$$

the average value of $1/r$ at $(0,0,1)$ is about .996, and the average at $(1.01,0,0)$ is about .991.

Python code below (please excuse the amateur-ness)

import random
import math

points = 10000000
e = .01
rad = 1+e
diam = 2*rad

pot= 0
count = 0
for i in range(1,points):
    x = diam*random.random()-rad
    y = diam*random.random()-rad
    z = diam*random.random()-rad
    r = math.sqrt((x-rad)*(x-rad)+y*y+z*z)
    if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
        pot = pot + 1/r
        count = count + 1

print pot/count

pot2 = 0
count = 0
for j in range(1,points):
    x = diam*random.random()-rad
    y = diam*random.random()-rad
    z = diam*random.random()-rad
    r = math.sqrt(x*x+y*y+(z-1.0)*(z-1.0))
    if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
        pot2 = pot2 + 1/r
        count = count + 1

print pot2/count

Answer 4

Start with the unperturbed gravitational potential for a uniform sphere of mass M and radius R, interior and exterior:

$$ \phi^0_\mathrm{in} = {-3M \over 2R} + {M\over 2R^3} (x^2 + y^2 + z^2) $$ $$ \phi^0_\mathrm{out} = {- M\over r} $$

Add a quadrupole perturbation, you get

$$ \phi_\mathrm{in} = \phi^0_\mathrm{in} + {\epsilon M\over R^3} D $$ $$ \phi_\mathrm{out} = \phi^0_\mathrm{out} + {M\epsilon R^2\over r^5} D $$

$$ D = x^2 + y^2 - 2 z^2 $$

The scale factors of M and R are just to make $\epsilon$ dimensionless, the falloff of $D\over r^5$ is just so that the exterior solution solves Laplace's equation, and the matching of the solutions is to ensure that on any ellipsoid near the sphere of radius R, the two solutions are equal to order $\epsilon$. The reason this works is because the $\phi^0$ solutions are matched both in value and in first derivative at x=R, so they stay matched in value to leading order even when perturbed away from a sphere. The order $\epsilon$ quadrupole terms are equal on the sphere, and therefore match to leading order.

The ellipsoid I will choose solves the equation:

$$ r^2 + \delta D = R^2 $$

The z-diameter is increased by a fraction $\delta$, while the x diameter decreased by $\delta/2$. So that the ratio of polar to equatorial radius is $3\delta/2$. To leading order

$$ r = R + {\delta D \over 2R}$$

We already matched the values of the inner and outer solutions, but we need to match the derivatives. taking the "d":

$$ d\phi_\mathrm{in} = {M\over R^3} (rdr) + {\epsilon M\over R^3} dD $$ $$ d\phi_\mathrm{out} = {M\over r^3} (rdr) + {MR^2\epsilon \over r^5} dD - {5\epsilon R^2 M\over r^7} (rdr) $$

$$ rdr = x dx + y dy + z dz $$ $$ dD = 2 x dx + 2ydy - 4z dz $$

To first order in $\epsilon$, only the first term of the second equation is modified by the fact that r is not constant on the ellipsoid. Specializing to the surface of the ellipsoid:

$$ d\phi_\mathrm{out}|_\mathrm{ellipsoid} = {M\over R^3} (rdr) + {3\delta \over 2 R^5}(rdr) + {\epsilon M \over R^3} dD - {5\epsilon M \over R^5} (rdr)$$

Equating the in and out derivatives, the parts proportional to $dD$ cancel (as they must--- the tangential derivatives are equal because the two functions are equal on the ellipsoid). The rest must cancel too, so

$$ {3\over 2} \delta = 5 \epsilon $$

So you find the relation between $\delta$ and $\epsilon$. The solution for $\phi_\mathrm{in}$ gives

$$ \phi_\mathrm{in} + {3M\over 2R} = {M\over 2R^3}( r^2 + {3\over 5} \delta D ) $$

Which means, looking at the equation in parentheses, that the equipotentials are 60% as squooshed as the ellipsoid.

Now there is a condition that this is balanced by rotation, meaning that the ellipsoid is an equipotential once you add the centrifugal potential:

$$ - {\omega^2\over 2} (x^2 + y^2) = -{\omega^2 \over 3} (x^2 + y^2 + z^2) -{\omega^2\over 6} (x^2 + y^2 - 2z^2) $$

To make the $\delta$ ellipsoid equipotential requires that $\omega^2\over 6$ equals the remaining ${2\over 5} {M\over 2R^2}$, so that, calling $M\over R^2$ (the acceleration of gravity) by the name "g", and $\omega^2 R$ by the name "C" (centrifugal)

$$\delta = {5\over 6} {C \over g} $$

The actual difference in equatorial and polar diameters is found by multiplying by 3/2 (see above):

$$ {3\over 2} \delta = {5\over 4} {C\over g} $$

instead of the naive estimate of ${C\over 2g}$. So the naive estimate is multiplied by two and a half for a uniform density rotating sphere.

Nonuniform interior: primitive model

The previous solution is both interior and exterior for a rotating uniform ellipsoid, and it is exact in r, it is only leading order in the deviation from spherical symmetry. So it immediately extends to give the shape of the Earth for a nonuniform interior mass distribution. The estimate with a uniform density is surprisingly good, and this is because there are competing effects largely cancelling out the correction for non-uniform density.

The two competing effects are: 1. the interior distribution is more elliptical than the surface, because the interior solution feels all the surrounding elliptical Earth deforming it, with extra density deforming it more. 2. The ellipticity of the interior is suppressed by the $1/r^3$ falloff of the quadrupole solution of Laplace's equation, which is $1/r^2$ faster than the usual potential. So although the interior is somewhat more deformed, the falloff more than compensates, and the effect of the interior extra density is to make the Earth more spherical, although not by much.

These competing effects are what shift the correction factor from 2.5 to 2, which is actually quite small considering that the interior of the Earth is extremely nonuniform, with the center more than three times as dense as the outer parts.

The exact solution is a little complicated, so I will start with a dopey model. This assumes that the Earth is a uniform ellipsoid of mass M and ellipticity parameter $\delta$, plus a point source in the middle (or a sphere, it doesn't matter), accounting for the extra mass in the interior, of mass M'. The interior potential is given by superposition. With the centrifugal potential:

$$ \phi_{int} = - {M'\over r} - {2M\over 3R} + {M\over 2R^3}(r^2 - {3\over 5} \delta D) + {\omega^2\over 2} r^2 - {\omega^2\over 6} D $$

This has the schematic form of spherical plus quadrupole (including the centrifugal force inside F and G)

$$ \phi_{int} = F(r) + G(r) D $$

The condition that the $\delta$ ellipsoid is an equipotential is found by replacing $r$ with $R - {\delta D\over 2R}$ inside F(r), and setting the D-part to zero:

$$ {F'(R) \delta \over 2R} = G(r) $$

In this case, you get the equation below, which reduces to the previous case when $M'=0$:

$$ {M'\over M+M'}\delta + {M\over M+M'} (\delta - {3\over 5} \delta) = - {C\over 3 g } $$

where $C=\omega^2 R$ is the centrifugal force, and $ g= {M+M'\over R^2} $ is the gravitational force at the surface. I should point out that the spherical part of the centrifugal potential ${\omega^2\over 2} r^2$ always contributes a subleading term proportional to $\omega^2\delta$ to the equation and should be dropped. The result is

$$ {3\over 2} \delta = {1\over 2 (1 - {3\over 5} {M\over M+M'}) } {C\over g} $$

So that if you choose M' to be .2 M, you get the correct answer, so that the extra equatorial radius is twice the naive amount of ${C\over 2g}$.

This says that the potential at the surface of the Earth is only modified from the uniform ellipsoid estimate by adding a sphere with 20% of the total mass at the center. This is somewhat small, considering the nonuniform density in the interior contains about 25% of the mass of the Earth (the perturbing mass is twice the density at half the radius, so about 25% of the total). The slight difference is due to the ellipticity of the core.

Nonuniform mass density II: exact solution

The main thing neglected in the above is that the center is also nonspherical, and so adds to the nonspherical D part of the potential on the surface. This effect mostly counteracts the general tendency of extra mass at the center to make the surface more spherical, although imperfectly, so that there is a correction left over.

You can consider it as a superposition of uniform ellipsoids of mean radius s, with ellipticity parameter $\delta(s)$ for $0<s<R$ increasing as you go toward the center. Each is uniform on the interior, with mass density $|\rho'(s)|$ where $\rho(s)$ is the extra density of the Earth at distance s from the center, so that $\rho(R)=0$. These ellipsoids are superposed on top of a uniform density ellipsoid of density $\rho_0$ equal to the surface density of the Earth's crust:

I will consider $\rho(s)$ and $\rho_0$ known, so that I also know $|\rho'(s)|$, it's (negative) derivative with respect to s, which is the density of the ellipsoid you add at s, and I also know:

$$ M(r) = \int_0^r 4\pi \rho(s) s^2 ds $$

The quantity $M(s)$ is ${1\over 4\pi}$ times the additional mass in the interior, as compared to a uniform Earth at crust density. Note that $M(s)$ is not affected by the ellipsoidal shape to leading order, because all the nested ellipsoids are quadrupole perturbations, and so contain the same volume as spheres.

Each of these concentric ellipsoids is itself an equipotential surface for the centrifugal potential plus the potential from the interior and exterior ellipsoids. So once you know the form of the potential of all these superposed ellipsoids, which is of the form of spherical + quadrupole + centrifugal quadrupole (the centrifugal spherical part always gives a subleading correction, so I omit it):

$$ \phi_\mathrm{int}(r) = F(r) + G(r) D + {\omega^2 \over 6} D $$

You know that each of these nested ellipsoids is an equipotential

$$ F(s - {\delta(s) \over 2s}) D + G(s) D - {\omega^2\over 6} D $$

so that the equation demanding that this is an equipotential at any s is

$$ {\delta(s) F'(s) \over 2s} - G(s) + {\omega^2\over 6} = 0 $$

To find the form of F and G, you first express the interior/exterior solution for a uniform ellipsoid in terms of the density $\rho$ and the radius R:

$$ {\phi_\mathrm{int}\over 4\pi} = - {\rho R^2\over 2} + {\rho\over 6} r^2 + {\rho \delta\over 10} D $$

$$ {\phi_\mathrm{ext}\over 4\pi} = - {\rho R^3 \over 3 r} + {\rho\delta R^5\over 10 r^5} D $$

You can check the sign and numerical value of the coefficients using the 3/5 rule for the interior equipotential ellipsoids, the separate matching of the spherical and D perturbations at r=R, and dimensional analysis. I put a factor $4\pi$ on the bottom of $\phi$ so that the right hand side solves the constant free form of Laplace's equation.

Now you can superpose all the ellipsoids, by setting $\delta$ on each ellipsoid to be $\delta(s)$, setting $\rho$ on each ellipsoid to be $|\rho'(s)|$, and $R$ to be $s$. I am only going to give the interior solution at r (doing integration by parts on the spherical part, where you know the answer is going to turn out to be, and throwing away some additive constant C) is:

$$ {\phi_\mathrm{int}(r)\over 4\pi} - C = {\rho_0\over 6} r^2 + {\rho_0 \delta(R)\over 10} D - {M(r)\over 4\pi r} + {1\over 10r^5} \int_0^r |\rho'(s)| \delta(s) s^5 ds D + {1\over 10} \int_r^R |\rho'(s)|\delta(s) D $$

The first two terms are the interior solution for constant density $\rho_0$. The third term is the total spherical contribution, which is just as in the spherical symmetric case. The fourth term is the the superposed exterior potential from the ellipsoids inside r, and the last term is the superposed interior potential from the ellipsoids outside r.

From this you can read off the spherical and quadrupole parts: $$ F(r) = {\rho_0\over 6} r^2 + {M(r)\over r} $$ $$ G(r) = {\rho_0\delta(R)\over 10} + {1\over 10r^5} \int_0^r |\rho'(s) |\delta(s) s^5 ds + {1\over 10} \int_r^R |\rho'(s)|\delta(s) $$

So that the integral equation for $\delta(s)$ asserts that the $\delta(r)$ shape is an equipotential at any depth.

$$ {F'(r)\delta(r)\over 2r} - G(r) + {\omega^2 \over 6} = 0 $$

This equation can be solved numerically for any mass profile in the interior, to find the $\delta(R)$. This is difficult to do by hand, but you can get qualitative insight.

Consider an ellipsoidal perturbation inside a uniform density ellipsoid. If you let this mass settle along an equipotential, it will settle to the same ellipsoidal shape as the surface, because the interior solution for the uniform ellipsoid is quadratic, and so has exact nested ellipsoids of the same shape as equipotentials. But this extra density will contribute less than it's share of elliptical potential to the surface, diminishing as the third power of the ratio of the radius of the Earth to the radius of the perturbation. But it will produce stronger ellipses inside, so that the interior is always more elliptical than the surface.

Oblate Core Model

The exact solution is too difficult for paper and pencil calculations, but looking [here]( http://www.google.com/imgres?hl=en&client=ubuntu&hs=dhf&sa=X&channel=fs&tbm=isch&prmd=imvns&tbnid=hjMCgNhAjHnRiM:&imgrefurl=http://www.springerimages.com/Images/Geosciences/1-10.1007_978-90-481-8702-7_100-1&docid=ijMBfCAOC1GhEM&imgurl=http://img.springerimages.com/Images/SpringerBooks/BSE%253D5898/BOK%253D978-90-481-8702-7/PRT%253D5/MediaObjects/WATER_978-90-481-8702-7_5_Part_Fig1-100_HTML.jpg&w=300&h=228&ei=ZccgUJCTK8iH6QHEuoHICQ&zoom=1&iact=hc&vpx=210&vpy=153&dur=4872&hovh=182&hovw=240&tx=134&ty=82&sig=108672344460589538944&page=1&tbnh=129&tbnw=170&start=0&ndsp=8&ved=1t:429,r:1,s:0,i:79&biw=729&bih=483 ), you see that it is sensible to model the Earth as two concentric spheres of radius $R$ and $R_1$ with total mass $M$ and $M_1$ and $\delta$ and $\delta_1$.

I will take

$$ R_1 = {R\over 2} $$

and

$$ M_1 = {M\over 4} $$

that is, the inner sphere is 3000 km across, with twice the density, which is roughly accurate. Superposing the potentials and finding the equation for the $\delta$s (the two point truncation of the integral equation), you find

$$ -\delta + {3\over 5} {M_0\over M_0 + M_1} \delta + {3\over 5} {M_1\over M_0 + M_1} \delta_1 ({R_1\over R})^2 = {C\over 3g} $$

$$ {M_0 \over M_0 + M_1} (-\delta_1 + {3\over 5} \delta) + {M_1 \over M_0 + M_1}( -\delta_1 + {3\over 5} \delta_1) = {C\over 3g} $$

Where

$$ g = {M_0+ M_1\over R^2}$$ $$ C = \omega^2 R $$

are the gravitational force and the centrifugal force per unit mass, as usual. Using the parameters, and defining $\epsilon = {3\delta\over 2}$ and $\epsilon_1={3\delta_1\over 2}$, one finds:

$$ - 1.04 \epsilon + .06 \epsilon = {C\over g} $$ $$ - 1.76 \epsilon_1 + .96 \epsilon = {C\over g} $$

(these are exact decimal fractions, there are denominators of 100 and 25). Subtracting the two equations gives:

$$ \epsilon_1 = {\epsilon\over .91} $$

(still exact fractions) Which gives the equation

$$ (-1.04 + {.06\over .91} ) \epsilon = {C\over g}$$

So that the factor in front is $.974$, instead of the naive 2. This gives an equatorial diameter of 44.3 km, as opposed to 42.73, which is close enough that the model essentially explains everything you wanted to know.

The value of $\epsilon_1$ is also interesting, it tells you that the Earth's core is 9% more eccentric than the outer ellipsoid of the Earth itself. Given that the accuracy of the model is at the 3% level, this should be very accurate.

Answer 5

In this answer, I will present a framework to use, and then I will frame the prior answers within that framework. Let me sum up the values we have here. I'll use the same notation (as best as possible) as everyone else and Wikipedia for an oblate spheroid where $a$ is the large, equatorial, radius.

• Mark1, method in the question, $2 (a-b) = 21.6 km$
• Mark2, method in the answer, past answer times 5/2 for $2 (a-b) = 54.0 km$

Here is my approach:

The mass of the earth can be taken a combination of two shapes, an inner sphere with radius $b$ and an edge volume, which is the oblate spheroid minus the inner sphere. The Earth has a certain average density defined by $\rho = M/V$, but this can be divided up into two different type of materials, the core, and the crust. The total mass requirement will then require that $M = V_{core} \rho_{core} + V_{edge} \rho_{edge}$, while the average density of the Earth requires $V = V_{core} + V_{edge}$, constraining $a$ and $b$ by one degree of freedom. When writing a code we can say that $a$ always implies $b$ according to $b=3V/(4 \pi a^2)$, the crust density also can be taken to imply the core density. Then the equipotential surface constraint dockets another degree of freedom, which can be used to iteratively find the value of $a$. Illustration:

The potential from the inner sphere is easy. I'll write it for a point on the equator and on the pole, combined with the other potentials.

$$U_a = G \frac{M_{core}}{a} + U_{a,edge} + \frac{1}{2} \omega^2 a^2$$ $$U_b = G \frac{M_{core}}{b} + U_{b,edge}$$

Obviously, the hard part is calculating the potential from the ridiculously irregular edge shape. Before that, however, it's important to think of the physical implications of looking at the problem this way. To begin with, what are relevant densities relevant to Earth? Here is the average density, and then the density within about 0 to 200 km of the surface.

$$rho_{Earth} = 5520. \frac{kg}{m^3}$$ $$rho_{crust} = 3400. \frac{kg}{m^3}$$

When we actually solve the problem, we will specify the crust density, and then that will imply a density of the core sphere. Is this method accurate? No. The main thing it misses is that the center sphere is not a spherically symmetric distribution of matter. The densities will theoretically stratify according to constant potential lines. In other words, if there were a high-density core of the Earth, it would also be an oblate spheroid. Because of that reason, introducing the inner sphere does miss some detail, but this model could still be pretty good.

Implementing this is a little tricky, as others have pointed out, but focusing the calculations on the edge volume helps a great deal. I also used a Monte Carlo method, but subtracted out the core volume. That is to say, I shot points randomly into the edge volume as efficiently as reasonable. In order to do this I used a weighted method, and it seems to have worked out okay. With 5 million iterations I found random variations in the calculation of the potential at $a$ to have a standard deviation of about $15 m$ of gravitational potential equivalent and around $2 m$ of gravitational potential equivalent at $b$. The reason for the higher deviation at $a$ is because it has more mass close to it, and the sampling was unbiased in the yz-plane for calculations of the $a$ potential and unbiased in the xy-plane for calculations of the $b$ potential. Anyway, as I iterated on the $a$ value, I set a tolerance for $100 m$, because this should be considerably higher than the statistical variation and it's good enough for calculation of the bulge. To recap, this is my method:

For calculation of potential at $a$

• Sample two random values for the $y$ and $z$ values between $0$ and $a$
• If these two values lie outside of the ellipse of $\frac{y^2}{a^2}+\frac{z^2}{b^2}<1$ then sample two new values and try again (gives something like a 30% loss of efficiency)
• Sample an x-value between the surface of the inner sphere and the outer spheroid. If the (y,z) pair doesn't fall within the inner sphere, sample between $0$ and the outer spheroid surface.
• Tally the potential between the sampled (x,y,z) point and (a,0,0), using the M_edge mass
• Repeat between (-x,y,z) and (a,0,0)
• Tally the weight for this sample as the distance between the two surfaces times two.
• The potential at (a,0,0) is then the total tallied potential divided by the total tallied weight.
• Repeat a similar method to find the potential at (0,0,b)
• Numerically root find to satisfy the equipotential condition discussed above.

I did this, and for the different values of the crust density, I got the following.

• $\rho_{crust} = 0$, implying gravitational field is insensitive to flatness, gets $2(a-b)= 21.8 \pm 0.1 km$
• $\rho_{crust} = 3400. \frac{kg}{m^3}$, a reasonable value for the crust density, gets $2(a-b) = 34.7 \pm 0.1 km$
• $\rho_{crust} = 5520. \frac{kg}{m^3}$, a fully homogenous density of the Earth gets $2(a-b) = 54.0 \pm 0.1 km$

I thought these are good results because the first and last of them come close to the prior answers within numerical error and the reasonable value for crust density gets closer to the real value of $42 km$.

If anyone is interested, I can look into putting the code for this on github or something similar. Otherwise, it's a little longer than the others posted here, so for now I'll avoid cluttering this space.

Answer 6

I look at two models of a "fat earth":

1. a spherically symmetric interior with an aspherical surface layer in hydrostatic equilibrium. This analysis generalizes from the constant density assumed in other answers and thereby exhibits the sensitivity of the flattening to the surface density. I compare the result to those of various other answers.
2. To estimate the effect of interior asphericity, I re-look at the case analyzed by Ron Maimon of two concentric constant density oblate spheres, both in hydrostatic equilibrium. My calculation shows that the effect of the core oblateness on the surface flattening is small, compared to the impact of a dense core on the average density.

For both analyses, I use the following result: consider a thin shell of material with density $\rho_s$, mean radius $r_o$, and thickness

$$\delta r_s(\theta)=-\frac{2}{3}f r_o P_2(\cos \theta) $$

(These are polar coordinates with $P_2 = \frac{1}{2} ( 3\cos^2 \theta -1)$, the second Legendre polynomial.) Note that $\Delta r_s = \delta r_s(\pi/2) - \delta r_s(0) = f r_0$, so $f$ is the flattening. It's easy to show that the mean radius is indeed $r_o$ and the net mass of the shell is 0. (Think of a surface mass density $\rho_s \delta r_s$, which is negative where $\delta r_s<0$. Ignore for the moment the aphysicality of negative mass; in practice this shell will be superimposed on a sphere.)

Then, to first order in f, the field generated by this shell is

$$ r > r_o: \phi_e(r,\theta)= \frac{2}{5} f \frac{G \rho_s V(r_o)}{r} \left( \frac{r_o}{r} \right)^2 P_2(\cos\theta) $$

$$ r < r_o: \phi_e(r,\theta)= \frac{2}{5} f \frac{G \rho_s V(r_o)}{r_o} \left( \frac{r}{r_o} \right)^2 P_2(\cos\theta) $$

(Here $ V(r) = \frac{4 \pi}{3} r^3$ is the volume of a sphere with radius $r$.) I include an outline of the calculation at the end.)

---

1. spherically symmetric interior + surface layer

Let the average density of the interior be $\rho_a=M/V(r_o)$ and the density at the surface $\rho_s$. Then, with the aspherical surface layer described above, the potential on the surface is (to first order in f):

$$ \phi_s = - \frac{GM}{r_{o}}\left( 1+ \frac{2}{3} f P_2(\cos \theta) \right) + \frac{2}{5} f \frac{G V(r_{o}) \rho_s}{r_{o}} P_2(\cos \theta) \\ - \frac{1}{3} \omega^2 r_{o}^2 \left( 1 - P_2(\cos \theta) \right) $$

where the first term includes the change in the spherical potential (the "$mgh$" term) the second is the quadrupole potential, and the last is the pseudo-potential of the rotation (written in terms of Legendre polynomials). Roughly speaking, at the poles the surface is closer to the center, and hence deeper in the potential well, but the reduction in nearby mass partly counteracts this effect.

For hydrostatic equilibrium, this potential must be independent of $\theta$, so:

$$ -\frac{2}{3} \frac{M}{r_o} f + \frac{2}{5} \frac{\rho_s V(r_o)}{r_o} f = -\frac{1}{3} \omega^2 r_o^2$$

Solving: $$ \left( 1 - \frac{3}{5} \frac{\rho_s}{\rho_a} \right) f = \frac{1}{2} \frac{\omega^2 r_o^2}{GM/r_o} = \frac{1}{2} \frac{\omega^2 r_o}{g} = \frac{q}{2} $$

or $$ \Delta r_s = f r_o = \frac{1}{1-\frac{3}{5} \frac{\rho_s}{\rho_a}} \frac{q}{2} r_o$$

Numbers: I'm using $r_o=6370$ km, $\rho_a = 5520$ kg/m$^3$, and $q=3.438 \cdot 10^{-3}$

Cases.

1. limit of 0 surface density $(\rho_s/\rho_a \rightarrow 0)$.
   $\Delta r_s = (q/2)r_o = 10.95 \text{ km} $
   No surface density implies no quadrupole field, so this result is the same as that calculated in the question, and also the same as in Professor Morin's 2004 text, section 9.4, problem 8. (Professor Morin assigned this problem 3 out of 4 difficulty stars; perhaps it should have been 4 out of 4?)
2. $\rho_s=\rho_a$, which includes the case of constant density analyzed by Qmechanic and Ron Maimon.
   $\Delta r_s = (q/2)r_o /(1-3/5) = (5/2) (q/2) r_o = 27.37 \text{ km} $, in agreement.
3. $\rho_s=3400$ kg/m$^3$, as in AlanSE's numerical calculation.
   $\Delta r_s = (q/2)r_o /[1-(3/5)(3400/5520)] = 1.586 (q/2) r_o = 17.37 \text{ km} $, again in agreement.
4. $\rho_s/\rho_a = 4/5$, which corresponds to Ron Maimon's extra core mass.
   $\Delta r_s = (q/2)r_o /[1-(3/5)(4/5)] = (25/13) (q/2) r_o = 1.923 (q/2) r_o =21.06 \text{ km} $.

   ---

2. Two superimposed oblate spheres, each with constant density

Following Ron Maimon's analysis, I now add a core mass. Notation:

1. I'm going to call the smaller radius sphere #1 (like Ron Maimon), with nominal radius $r_{co}$, mass $M_1$, constant density $\rho_1=M_1/V(r_{co})$, and flattening $f_c$.
2. The larger radius sphere is #2 (different from Ron Maimon's analysis, sorry), with nominal radius $r_{so}$, mass $M_2$, constant density $\rho_2=M_2/V(r_{so})$, and flattening $f_s$

Both sphere surfaces are assumed to be in hydrostatic equilibrium (and hence oblate). The two flattening coefficients couple via their quadrupole fields.

Superposing the fields for the two masses, the potential at the core surface is:

$$ \phi_c = - \frac{G M_1}{r_{co}} \left( 1+ \frac{2}{3} f_c P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_1}{r_{co}} \left( \frac{r_{co}}{r_{co}} \right)^2 f_c P_2(\cos \theta) -\frac{3}{2} \frac{G M_2}{r_{so}} \left[ 1 - \frac{1}{3} \left( \frac{r_{co}}{r_{so}} \right)^2 \left( 1 - \frac{4}{3} f_c P_2(\cos \theta) \right) \right] + \frac{2}{5} \frac{G M_2}{r_{so}} \left( \frac{r_{co}}{r_{so}} \right)^2 f_s P_2(\cos \theta) - \frac{1}{3} \omega^2 r_{co}^2 \left( 1 - P_2(\cos \theta) \right) $$

The potential at the outer surface is:

$$ \phi_s = - \frac{G M_1}{r_{so}} \left( 1+ \frac{2}{3} f_s P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_1}{r_{so}} \left( \frac{r_{co}}{r_{so}} \right)^2 f_c \, P_2(\cos \theta) - \frac{G M_2}{r_{so}} \left( 1+ \frac{2}{3} f_s P_2(\cos \theta) \right) + \frac{2}{5} \frac{G M_2}{r_{so}} \left( \frac{r_{so}}{r_{so}} \right)^2 f_s \, P_2(\cos \theta) - \frac{1}{3} \omega^2 r_{so}^2 \left( 1 - P_2(\cos \theta) \right) $$

Requiring, for hydrostatic equilibrium, that there be no $\theta$-dependence gives the pair of equations:

$$ \left[ 1 - \frac{3}{5} \frac{M_2}{M_1+M_2} \right] f_s - \alpha f_c = \frac{q}{2} \text{ , surface} $$

$$ -\frac{3}{5} \frac{M_2}{M_1+M_2} f_s + \beta f_c = \frac{q}{2} \text{ , core} $$

$$ \text{where } \alpha = \frac{3}{5} \left( \frac{r_{co}}{r_{so}} \right)^2 \frac{M_1}{M_1+M_2} $$

$$ \text{and } \beta = \left( \frac{r_{so}}{r_{co}} \right)^3 \left[ \frac{M_1 + \left( \frac{r_{co}}{r_{so}} \right)^3 M_2}{M_1+M_2} - \frac{3}{5} \frac{M_1}{M_1+M_2} \right] $$

Solving: $$ f_c = \frac{1}{\beta + \alpha} f_s $$ $$ f_s = \frac{1}{1-\frac{3}{5} \frac{M_2}{M_1+M_2} - \frac{1}{1+\beta/\alpha}} \frac{q}{2} $$

Note that $\frac{M_2}{M_1+M_2} = \frac{\rho_s}{\rho_a}$, consistent with the previous analysis.

Numerics. Using Ron Maimon's values: $\frac{r_{co}}{r_{so}}=\frac{1}{2}, and \frac{M_1}{M_2}=\frac{1}{4}$, one calculates $\alpha = 0.030, \beta = 1.440 $ and $1/(1+\beta/\alpha) = 0.020$, with the result that the value for $f_s/(q/2)$ calculated in part 1 case 4, $1.923$, is increased to 2.002 (+4%) when the core oblateness is included. This effect is smaller than the effect of the core mass in lowering $\rho_s/\rho_a$, which reduced $f_s/(q/2)$ from 2.5 to 1.923 (-23%).

---

Finally, here is an outline of the calculation of the potential due to a "quadrupole" shell: approximate the shell as a surface of variable mass density at the mean radius $r_o$. This surface can be decomposed into rings of constant $\theta$ with mass

$$ dm(\theta') = 2 \pi r_0 \sin \theta' \, r_0 \, d\theta' \, \delta r_s \, \rho_s = - f \rho_s \frac{4 \pi}{3} r_o^3 \, P_2(\cos \theta') \, \sin \theta' \, d \theta'$$

The potential generated by this ring is calculated in Jackson, Classical EM, Section 3.3:

$$r>r_o: d\phi (r,\theta) = - \frac{G \, dm(\theta')}{r} \sum_{l=0}^{\infty} \left( \frac{r_o}{r} \right)^l P_l(\cos \theta') P_l(\cos \theta) $$

Substituting for the mass density (which is proprotional to $P_2$) and integrating over $\theta'$ zeros all terms except the $P_2$ term, due to the orthogonality of Legendre polynomials, giving the stated result. (The development for $r<r_0$ is similar.)

Answer 7

Here I would like to numerically check the theoretical prediction of a factor $\frac{2}{5}$ in difference from Mark Eichenlaub's original monopole argument. In practice, this means calculating the difference in the gravitational potential between the North pole and the Equator, and divide by the corresponding difference in monopole potentials. For numerical reasons, it is better in practice to calculate the inverse (=reciproc) fraction, which then should be compared with $\frac{5}{2}$. Since my programming skills are limited, I have just written a slow MAPLE code to do the job.

 b:=100; f:=.10; a:=b*(1+f); V1:=evalf(4*Pi*a^2*b/3); 
 xa:=a; ya:= 0; za:=0; xb:=0; yb:=0; zb:=b; 
 U1a := evalf(V1/sqrt(xa^2 + ya^2 + za^2)); 
 U1b := evalf(V1/sqrt(xb^2 + yb^2 + zb^2));

 Ua:=0;Ub:=0;V:=0;
 for x from -a-.5 by 1 to a+.5 do 
 for y from -a-.5 by 1 to a+.5 do 
 for z from -b-.5 by 1 to b+.5 do 
 if (x/a)^2 + (y/a)^2 + (z/b)^2 < 1 then 
 Ua:=Ua + 1/sqrt((x-xa)^2 + (y-ya)^2 + (z-za)^2); 
 Ub:=Ub + 1/sqrt((x-xb)^2 + (y-yb)^2 + (z-zb)^2); 
 V:=V+1;
 end if;od;od;od;

 b;f;Ua;U1a;Ub;U1b;V;V1;Ub-Ua;U1b-U1a;(U1b-U1a)/(Ub-Ua);

The result with polar radius $b:=100$ was

$$ \begin{array}{cc} \mathrm{Flatness}&\mathrm{Inverse \ fraction} \\ f & \frac{U_{1b}-U_{1a}}{U_b-U_a} \\ \\ 0.10 & 2.5549 \\ 0.05 & 2.5353 \\ 0.03 & 2.5149 \\ 0.02 & 2.5085 \\ 0.01 & 2.5088 \\ \end{array} $$

The fact that the estimate does not improve from flatness $0.02$ to $0.01$ is a lattice artifact, because the lattice spacing is of same order as the difference in radius between the North pole and the Equator.