The squared angular velocity needed for equilibrium reaches a maximum of 
In the interior of a solid ball of matter with uniform density, the Newtonian gravitational potential takes the simple form:
U(r) = 2 π G ρ (⅓ r2 – R2)
where G is the Newtonian gravitational constant, ρ is the density, and R is the radius of the ball.
Of course this is just a toy model that is quite different from real astronomical bodies, whose density will usually vary significantly with depth unless they are quite small — and if they are small they are less likely to be spherical.
What about an ellipsoid, rather than a spherical ball? It turns out that the potential in the interior of an ellipsoid with uniform density again takes on about the simplest form possible, given that the shape is no longer radially symmetrical:
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
We will prove this result and find the values of the four coefficients in this formula in terms of the dimensions of the ellipsoid. We will then examine some of the consequences of this for slightly more complicated systems, such as a uniformly spinning body made of a fluid of constant density.
Most of the results discussed here are centuries old, and I will not make any attempt to give detailed attributions. There is a wonderful survey of the whole topic by Chandrasekhar[1] that traces some of the history as well as delving deeply into the mathematical physics.
Consider a shell of uniform density, whose boundaries are two concentric spheres. Newton gave a simple argument that a body located anywhere inside this shell will experience no gravitational force from the shell.
Suppose the body is located at point P in the diagram on the left. If we consider the material of the shell that lies within a narrow double cone with its apex at P, it can be subdivided into slices that lie at an equal distance r from P. The amount of matter in each such slice will be proportional to r2 dr, but the inverse-square law of gravity means the force at P is proportional to 1/r2, and hence the total force from each slice is independent of r and just proportional to its thickness, dr.
This means the total force at P in the direction of B is proportional to the distance AB, while the total force in the direction of C is proportional to the distance CD. But because the boundaries of the shell are concentric spheres we always have AB = CD, so there are equal and opposite forces at P that (in the limit of an infinitesimal cone) cancel out exactly. By summing over infinitesimal cones and infinitesimal slices, we can account for all the matter in the shell this way, so the effect of the entire shell is that there is zero net force, for any choice of P.
If there is no gravitational force inside the shell, the gravitatitional potential, U, must be constant. That means the potential everywhere inside the shell will be equal to the potential at the centre.
For an infinitesimally thin shell with a radius S and a uniform surface density of σ, the potential at the centre, and hence the potential throughout the interior, is:
Ushell, int(S, σ, r)
= –G M / S
= –G (4 π S2 σ) / S
= –4 π G σ S
We can confirm (and extend) this result with a direct computation. The potential at some point P at a distance r from the centre of an infinitesimal shell of radius S and surface density σ can be found by integrating over the surface.
Ushell(S, σ, r) = –2 π G σ S2 ∫0π sin θ / √[S2 + r2 – 2 S r cos θ] dθ
Here θ is the angle between a ray from the centre of the shell that passes through the point P and a ring of matter on the shell that is equidistant from P. The ring of matter has a radius of S sin θ and a thickness measured along the surface of the shell of S dθ, and the distance from P is given by the denominator in the integral above.
With a change of variable to u = cos θ we have:
Ushell(S, σ, r) = –2 π G σ S2 ∫–11 1 / √[S2 + r2 – 2 S r u] du
= (2 π G σ S / r) (√[S2 + r2 – 2 S r u]) |–11
= (2 π G σ S / r) (√[S2 + r2 – 2 S r] – √[S2 + r2 + 2 S r])
= (2 π G σ S / r) (|S – r| – |S + r|)
This gives us, for the interior of the shell where r ≤ S, and the exterior where r ≥ S:
Potential for an infinitesimal spherical shell of radius S and uniform surface density σ
Ushell, int(S, σ, r) = –4 π G σ S
Ushell, ext(S, σ, r) = –4 π G σ S2 / r
So, as we showed previously, the potential inside the shell is constant (and hence a test particle there would experience no gravitational force) while the potential outside the shell is the same as that due to a point mass (of the same mass as the whole shell) located at the centre of the shell.
We can now compute the potential inside or outside a solid ball of uniform density ρ and radius R, by integrating the contribution from a continuum of shells whose radius S ranges from 0 to R.
If the point of interest lies at a distance r from the centre of the ball, then for the interior case, where r ≤ R, the point lies in the exterior of all the shells with radius ranging from 0 to r, but in the interior of those with radius ranging from r to R:
Uball, int(R, ρ, r) = ∫0r Ushell, ext(S, ρ, r) dS + ∫rR Ushell, int(S, ρ, r) dS
= –4 π G ρ (∫0r (S2 / r) dS + ∫rR S dS)
= –4 π G ρ (⅓ r2 + ½ R2 – ½ r2)
= 2 π G ρ (⅓ r2 – R2)
For the exterior case, where r ≥ R, the point of interest lies in the exterior of all the shells:
Uball, ext(R, ρ, r) = ∫0R Ushell, ext(S, ρ, r) dS
= –4 π G ρ ∫0R (S2 / r) dS
= –(4/3) π G ρ R3 / r
In the exterior case, the potential is the same as that of a point mass, of the same total mass as the ball, located at its centre.
In the interior case, the potential increases with the square of the distance from the centre, so the gravitational acceleration, which is the opposite of the gradient of the potential, is proportional to the distance from the centre. This means that a test particle dropped into a borehole in a solid sphere of uniform density would undergo simple harmonic motion like a weight on the end of a spring, a famous thought experiment that you can read much more about here.
Potential for a solid ball of radius R and uniform density ρ
Uball, int(R, ρ, r) = 2 π G ρ (⅓ r2 – R2)
Uball, ext(R, ρ, r) = –(4/3) π G ρ R3 / r
The same argument that shows that there is no net gravitational force on a body inside a spherical shell of uniform density works just as well for an ellipsoidal shell of uniform density — so long as the inner and outer boundaries of the shell are concentric, aligned along the same axes, and are similar ellipsoids (that is, they are the same shape apart from an overall scale factor).
If the shell meets those requirements, we can be sure that the distance AB in the diagram on the left will always equal the distance CD, because we can obtain the ellipsoidal shell by applying a suitable linear transformation to a spherical shell. A linear transformation will always multiply the lengths of any two parallel line segments by the same factor, so the pairs of equal distances from the spherical case will remain equal.
A shell like this is known as a homoeoid. Note that the thickness of a homoeoid, measured orthogonally to the surface, is not constant. This means that an infinitesimally thin ellipsoidal shell with a constant surface density, i.e. a fixed quantity of mass per unit area, will not be a homoeoid, and it will not exert zero net gravitational force in its interior. Of course we can still study infinitesimally thin homoeoids, but they will not have constant surface densities.
An ellipsoid with semi-axes of length a, b and c has a boundary that can be described (in a coordinate system with its origin at the centre of the ellipsoid and its axes aligned with those of the ellipsoid) by the quadratric equation:
x2/a2 + y2/b2 + z2/c2 = 1
We will assume that our ellipsoid has a uniform density of ρ.
As with the ellipsoidal shell, we will consider the gravitational force at some point P due to the matter in a narrow double cone with its apex at P, but we will consider two such double cones at the same time, with the y- and z-components of the two cones’ direction vectors being the opposite of each other, as in the diagram on the left.
As we noted before, the amount of material in each slice of the cone at a given distance r is proportional to r2 dr, while the inverse-square law means the force is proportional to 1/r2. So the force due to each slice is just proportional to its thickness, dr, and the total force at P due to the matter in each cone will be proportional to the length of the cone.
It follows that the x-component of the force at P due to each cone is proportional to the extent of the cone in the x direction, Δx, where this is a signed quantity, with a force in the negative x direction when Δx is negative.
If P has coordinates (x, y, z), the four values of Δx are the roots of the two quadratic equations:
Q±(Δx) = (x + Δx)2 / a2 + (y ± β Δx)2 / b2 + (z ± γ Δx)2 / c2 – 1 = 0
The parameters β and γ determine the direction of the cones, while the sign determines which of the two double cones we are looking at.
We don’t need to explicitly solve these quadratics, because we don’t care about the individual values of Δx, just the sum of all four values. The sum of the two roots of any quadratic is simpler than the individual roots:
If Q2 Δx2 + Q1 Δx + Q0 = 0
then Δx1 + Δx2 = –Q1/Q2
For our particular quadratics, we have:
Q2 = 1 / a2 + β2 / b2 + γ2 / c2
Q1,± = 2 x / a2 ± 2 β y / b2 ± 2 γ z / c2
Δx1,+ + Δx2,+ + Δx1,– + Δx2,–
= –(Q1,+ + Q1,–) / Q2
= –4 x / (a2 Q2)
So the total x-component of the force at P due to the matter in all four cones is independent of y and z, and directly proportional to x. The constant of proportionality will vary with the direction of the cones, but it will not depend on P at all, so we can integrate over all directions for the cones, taking in all the material of the ellipsoid, to obtain a single constant.
The same analysis can be applied to the y-component and z-component of the force, with a different constant of proportionality for each axis. We will compute the constants themselves later, by a slightly different method, but what we have shown is that the gravitational force at P has components along each axis that are proportional to the displacement from the centre along that axis.
From that simple, linear relationship, it follows that the gravitational potential must take a simple quadratic form, which we can write as:
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
Instead of trying to compute the gravitational potential at an arbitrary point in the interior directly, we will perform a couple of simpler calculations: finding the potential at the centre of the ellipsoid, which we will call U0, and the potentials at the vertices of the ellipsoid, where the axes intersect the boundary, which we will call Ua, Ub and Uc. Those four quantities will be all we need to find the four coefficients A, B, C and D.
To compute U0, imagine sweeping a vector s = (sx, sy, sz) across a unit sphere with the same centre as the ellipsoid, and summing the potential due to the matter in a cone around s that subtends an infinitesimal solid angle ds and reaches from the centre to the boundary. The height of the cone, h0(s), will satisfy:
h0(s)2 sx2/a2 + h0(s)2 sy2/b2 + h0(s)2 sz2/c2 = 1
h0(s)2 = 1 / (sx2/a2 + sy2/b2 + sz2/c2)
The potential at the centre due to the matter in each cone can be found by integrating over the distance from the centre, r, from 0 to h0(s):
dU0 = (–G ρ ∫0h0(s) r2/r dr) ds
= –½ G ρ h0(s)2 ds
If we adopt standard polar coordinates where:
s = (sin θ cos φ, sin θ sin φ, cos θ)
ds = sin θ dφ dθ
we then have an integral for U0:
U0(a, b, c) = –½ G ρ ∫0π ∫02π sin θ / [sin2 θ cos2 φ / a2 + sin2 θ sin2 φ / b2 + cos2 θ / c2] dφ dθ
Before computing this for the generic case where none of the semi-axes are equal, suppose we have an ellipsoid of revolution, with two equal semi-axes, and choose the axes so that a = b. The integrand then becomes independent of φ, so we can replace that part of the integration by a factor of 2π, giving us:
U0(a, a, c) = –π G ρ ∫0π sin θ / [sin2 θ / a2 + cos2 θ / c2] dθ
If a > c (an oblate spheroid), we can make the substitution:
tan α = √[a2/c2 – 1] cos θ
sec2 α dα = (1 + tan2 α) dα = –√[a2/c2 – 1] sin θ dθ
(1 + (a2/c2 – 1) cos2 θ) dα = –√[a2/c2 – 1] sin θ dθ
a2 [sin2 θ / a2 + cos2 θ / c2] dα = –√[a2/c2 – 1] sin θ dθ
So we can integrate with respect to α easily and obtain:
U0(a, a, c) = (–2π G ρ a2 / √[a2/c2 – 1]) atan(√[a2/c2 – 1])
The atan here can be expressed a bit more simply as an acos, and a similar substitution with tanh in place of tan and a change of sign inside the square root does the same job for a prolate spheroid.
Potential at the centre of an ellipsoid of revolution with uniform density ρ, semi-axes a, a, c
For a > c, oblate spheroid
U0(a, a, c) = (–2π G ρ a2 c / √[a2 – c2]) acos(c/a)
For a < c, prolate spheroid
U0(a, a, c) = (–2π G ρ a2 c / √[c2 – a2]) acosh(c/a)
In fact either formula is mathematically correct for both the oblate and prolate cases, if we allow intermediate values to be imaginary numbers with factors of i that cancel out, but for the purpose of evaluating the result in programming languages where it’s easier to work entirely with real numbers it’s worth having these two distinct formulas.
In the generic case we need to perform an integral over the polar coordinate φ, of the form:
I(p, q) = ∫0π/2 dφ / (p + q cos2 φ)
This can be evaluated with the substitution:
tan φ = √[1 + q/p] tan γ
sec2 φ = 1 + tan2 φ = 1 + (1 + q/p) tan2 γ
p sec2 φ + q = (p + q) (1 + tan2 γ) = (p + q) sec2 γ
sec2 φ dφ = √[1 + q/p] sec2 γ dγ
dφ / (p + q cos2 φ)
= √[1 + q/p] sec2 γ / [sec2 φ (p + q cos2 φ)] dγ
= √[1 + q/p] sec2 γ / [p sec2 φ + q] dγ
= √[1 + q/p] / (p + q) dγ
I(p, q) = π/2 / √[p (p + q)]
Note that if the range of integration for φ is extended from π/2 to 2π we just get four times the original integral.
We can put the φ part of the generic integral for
p = cos2 θ / c2 + sin2 θ / b2
q = sin2 θ (1/a2 – 1/b2)
Using our result for I(p, q), applying sin2 θ = 1 – cos2 θ, and moving some factors outside the integral, we get:
U0(a, b, c) = –π G ρ a b ∫0π sin θ / √[(1 – (1–b2/c2) cos2 θ) (1 – (1–a2/c2) cos2 θ)] dθ
With the substitution w = cos θ and using the symmetry with respect to the sign of w this becomes:
U0(a, b, c) = –2 π G ρ a b ∫01 1 / √[(1 – (1–b2/c2) w2) (1 – (1–a2/c2) w2)] dw
We can rewrite this in the form of an incomplete elliptic integral of the first kind, a class of functions defined by:
F(φ | m) = ∫0sin(φ) dt / √[(1–t2) (1–mt2)]
where –π/2 < φ < π/2 and m sin2 φ < 1
Putting w = t / √[1–b2/c2] and using asin(√[1–b2/c2]) = acos(b/c), we arrive at:
Potential at the centre of an ellipsoid with uniform density ρ, semi-axes a < b < c
U0(a, b, c) = (–2 π G ρ a b c / √[c2–b2]) F(acos(b/c) | (c2–a2) / (c2–b2))
If we choose a < b < c, then the formula above will involve only real numbers. Other permutations of a, b and c can result in imaginary numbers appearing in intermediate steps in the calculation. So, in all our formulas for generic ellipsoids we will adopt a strict ordering of a < b < c.
There is an alternative form for U0 that is manifestly symmetrical in the semi-axes. It is not as useful for calculations, as the integral remains unevaluated, but it will make some mathematical properties easier to demonstrate. If we go back to the integral over θ, split it into two equal parts with endpoints at 0 and π/2, and make the substitution u = c2 tan2 θ, we obtain:
Potential at the centre of an ellipsoid with uniform density ρ, symmetrical form
U0(a, b, c) = (–π G ρ a b c) ∫0∞ du / Δ
where Δ = √[(a2 + u)(b2 + u)(c2 + u)]
To compute the potential at a vertex of the ellipsoid, such as the point (0, 0, a), we shift the apex of the cones that sweep out the volume of the ellipsoid from the centre to the vertex, and the height of the cones, ha(s), satisfies a different equation:
(a – ha(s) sx)2/a2 + ha(s) sy2/b2 + ha(s) sz2/c2 = 1
ha(s) = (2 sx / a) / (sx2/a2 + sy2/b2 + sz2/c2)
We then have:
dUa = –½ G ρ ha(s)2 ds
but we only need to sweep s over half of the unit sphere, because half the directions at the vertex just hit the boundary immediately.
However, it turns out that we do not need to perform this integration at all! If we take the derivative of the integrand we used to evaluate U0 with respect to one of the semi-axis lengths, we find:
∂a (h0(s)2)
= ∂a (1 / (sx2/a2 + sy2/b2 + sz2/c2))
= (2 sx2/a3) / (sx2/a2 + sy2/b2 + sz2/c2)2
= ha(s)2 / (2a)
We can use this to compute Ua by swapping the order of integration and differentiation. Taking account of the fact that for the vertex potential we only integrate over the half the unit sphere pointing one way along one of the axes, which for the potential at the centre simply halves the result, this gives us:
Ua = a ∂aU0
and analogous formulas for Ub and Uc. The only complication is that for the ellipsoids of rotation, if we want the potential at the vertex for one of the semi-axes whose length is identical to another, we need to halve the formula given here, because the partial derivative applied to a single variable that gives the lengths of both semi-axes varies them both, and yields double the result.
Potential at the vertices of an ellipsoid of revolution with uniform density ρ, semi-axes a, a, c
For a > c, oblate spheroid
Ua(a, a, c) = (–π G ρ a2 c) ((a2 – 2c2) acos(c/a) / (a2 – c2)3/2 + c / (a2 – c2))
Uc(a, a, c) = (–2π G ρ a2 c) (a2 acos(c/a) / (a2 – c2)3/2 – c / (a2 – c2))
For a < c, prolate spheroid
Ua(a, a, c) = (–π G ρ a2 c) ((2c2 – a2) acosh(c/a) / (c2 – a2)3/2 – c / (c2 – a2))
Uc(a, a, c) = (–2π G ρ a2 c) (–a2 acosh(c/a) / (c2 – a2)3/2 + c / (c2 – a2))
Potential at the vertices of an ellipsoid with uniform density ρ, semi-axes a < b < c
Ua(a, b, c) = –2π G ρ a b c (–a b / (c (b2–a2)) + a2 E √[c2–b2] / ((b2–a2) (c2–a2)) + c2 F / ((c2–a2) √[c2–b2]))
Ub(a, b, c) = –2π G ρ a b c (a b / (c (b2–a2)) – b2 E / ((b2–a2) √[c2–b2]) + F/√[c2–b2])
Uc(a, b, c) = –2π G ρ a b c (c2 E – a2 F) / ((c2–a2) √[c2–b2])
Here F = F(acos(b/c) | (c2–a2) / (c2–b2)) and E = E(acos(b/c) | (c2–a2) / (c2–b2))
are incomplete elliptic integrals of the first and second kind.
Potential at the vertices of an ellipsoid with uniform density ρ, symmetrical form
Ua(a, b, c) = U0(a, b, c) + (π G ρ a3 b c) ∫0∞ du / ((a2 + u) Δ)
Ub(a, b, c) = U0(a, b, c) + (π G ρ a b3 c) ∫0∞ du / ((b2 + u) Δ)
Uc(a, b, c) = U0(a, b, c) + (π G ρ a b c3) ∫0∞ du / ((c2 + u) Δ)
where Δ = √[(a2 + u)(b2 + u)(c2 + u)]
Now suppose we want to find the potential at some arbitrary point P = (x, y, z) in the interior of the ellipsoid. The potentials we have already computed, at the centre of the ellipsoid and at the three vertices, are enough to determine the values of the four coefficients A, B, C and D.
A = (Ua – U0) / (2 π G ρ a2)
B = (Ub – U0) / (2 π G ρ b2)
C = (Uc – U0) / (2 π G ρ c2)
D = U0 / (2 π G ρ)
Making use of our results for U0, Ua, Ub and Uc, we obtain the complete solutions:
Potential in the interior of an ellipsoid of revolution with uniform density ρ, semi-axes a, a, c
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
For a > c, oblate spheroid
A = B = –½c2 / (a2–c2) + ½a2 c acos(c/a) / (a2–c2)3/2
C = a2 / (a2–c2) – a2 c acos(c/a) / (a2–c2)3/2
D = –a2 c acos(c/a) / √[a2–c2]
For a < c, prolate spheroid
A = B = ½c2 / (c2–a2) – ½a2 c acosh(c/a) / (c2–a2)3/2
C = –a2 / (c2–a2) + a2 c acosh(c/a) / (c2–a2)3/2
D = –a2 c acosh(c/a) / √[c2–a2]
Potential in the interior of an ellipsoid with uniform density ρ, semi-axes a < b < c
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
A = b2/(b2–a2) – a b c (E √[c2–b2] / ((b2–a2) (c2–a2)) + F / ((c2–a2) √[c2–b2]))
B = –a2/(b2–a2) + a b c E / ((b2–a2) √[c2–b2])
C = a b c (F – E) / ((c2–a2) √[c2–b2])
D = –a b c F/√[c2–b2]
Here F = F(acos(b/c) | (c2–a2) / (c2–b2)) and E = E(acos(b/c) | (c2–a2) / (c2–b2))
are incomplete elliptic integrals of the first and second kind.
Potential in the interior of an ellipsoid with uniform density ρ, symmetrical form
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
A = ½ a b c ∫0∞ du / ((a2 + u) Δ)
B = ½ a b c ∫0∞ du / ((b2 + u) Δ)
C = ½ a b c ∫0∞ du / ((c2 + u) Δ)
D = –½ a b c ∫0∞ du / Δ
where Δ = √[(a2 + u)(b2 + u)(c2 + u)]
Acceleration due to gravity in the interior of an ellipsoid with uniform density
g(x, y, z) = –∇U(x, y, z)
= –4 π G ρ (A x, B y, C z)
We’ve derived all our potentials by integrating dU = –G ρ dV / r over the matter responsible for the gravitational field, but there is also a local differential equation that the potential must satisfy. Poisson’s equation for gravity states that the Laplacian operator applied to the Newtonian gravitational potential is equal to 4π G ρ. Applying this to our quadratic for U, we get:
U(x, y, z) = 2 π G ρ (A x2 + B y2 + C z2 + D)
∂2U/∂x2 + ∂2U/∂y2 + ∂2U/∂z2 = 4π G ρ (A + B + C)
So as a check on our results, we note that for all our sets of coefficients we have, as required:
A + B + C = 1
The symmetrical form of the expressions for the coefficients A, B and C lets us see fairly easily how some quantities will differ when they are associated with either of two unequal semi-axes. The integrands for any two of these coefficients are always positive, and the ratio between them is:
Integrand ratio A on B = (b2 + u) / (a2 + u)
If two integrands are in a ratio that is always less than 1, or always greater than 1, for all nonnegative values of u (except perhaps for an isolated point), then the same will be true of the results of the integration. So if a < b, we can easily show that A > B, and that A a2 < B b2. From the latter, it follows that Ua < Ub, i.e. the vertex closer to the centre of the ellipsoid is “deeper” in the potential well.
It’s a bit trickier to determine the ordering of A a and B b, which tells us the relative strength of the gravitational acceleration at the vertices. The two integrands swap their relative size at u = a b, and it is not obvious which part of the integral dominates.
Inequalities that hold if a < b
Integrand ratio A on B = (b2 + u) / (a2 + u) > 1 for all u ≥ 0
So A > B
Integrand ratio (A a2) on (B b2) = (a2 b2 + a2 u) / (a2 b2 + b2 u) < 1 for all u > 0
So A a2 < B b2 and Ua < Ub
We have computed the potential and the acceleration due to gravity in the interior of an ellipsoid with uniform density, but what can we say about the gravitational field in the surrounding vacuum? The key to this turns out to be a simple result about the geometry of confocal ellipsoids.
What are confocal ellipsoids? Let’s start with the idea of confocal ellipses, in two dimensions. Every ellipse has two points known as its foci, which lie on the major axis on either side of the centre. In the standard orientation where the major axis lies on the x-axis and the semi-axes of the ellipse are a and b, the foci are located at (±f, 0), where f 2 = a2 – b2.
It follows that if we change a and b in such a way that a2 and b2 are both changed by the same amount, say λ, then f will be unchanged, and the foci will be the same. So we have a family of confocal ellipses, described by the family of equations where we fix a and b and allow λ to take on different values:
x2 / (a2 + λ) + y2 / (b2 + λ) = 1
Similarly, we have a family of confocal ellipsoids, described by the family of equations for fixed a, b, c and different choices of λ:
x2 / (a2 + λ) + y2 / (b2 + λ) + z2 / (c2 + λ) = 1
Given two ellipsoids, E1 and E2, with semi-axes ax, ay, az and αx, αy, αz, we can map points from the surface of E1 to the surface of E2 by rescaling each coordinate as follows:
T(x, y, z) = ((αx /ax) x, (αy /ay) y, (αz /az) z)
If E1 and E2 are confocal, then for μ = x, y, z we have:
(αμ/aμ)2 – 1
= (αμ2 – aμ2) / aμ2
= λ / aμ2
for some fixed λ that does not depend on the particular axis.
Let P = (px, py, pz) and Q = (qx, qy, qz) be any two points on the surface of E1.
dist(P, T(Q))2 – dist(T(P), Q)2
= Σμ [(pμ – (αμ/aμ) qμ)2 – ((αμ/aμ) pμ – qμ)2]
= Σμ [(qμ2 – pμ2) ((αμ/aμ)2 – 1)]
= λ Σμ [(qμ2 – pμ2) / aμ2]
= λ (Σμ [qμ2 / aμ2] – Σμ [pμ2 / aμ2])
= λ (1 – 1)
= 0
So dist(P, T(Q)) = dist(T(P), Q).
Another result we will need is that for a rod of uniform density, the component of the gravitational attraction of the rod, along the direction of the rod itself, at some point P, is proportional to the difference of the reciprocals of the distances of P from the endpoints of the rod, A and B. If we suppose the rod lies along the x-axis with endpoints at xA and xB and has cross-sectional area σ, while P lies at (0, d, 0), we have:
gx = G ρ σ ∫xAxB x / (d2 + x2)3/2 dx
= G ρ σ [1/(d2 + xA2)½ – 1/(d2 + xB2)½]
= G ρ σ [1/dist(P, A) – 1/dist(P, B)]
Now, consider two confocal ellipsoids, E1 and E2, with points P on the first and T(P) on the second, where T is the map between ellipsoids we defined previously. Let E1 be the smaller of the two ellipsoids. First, suppose that E1 is a solid ellipsoid of uniform density ρ, and we compute the gravitational attraction g(E1, T(P)) it exerts at the exterior point T(P). The point T(P) lies on E2, but here E2 is just a geometric construct to aid our calculations, containing no material. Then suppose we swap the roles of the ellipsoids, with E2 the solid ellipsoid of uniform density ρ, and we compute the gravitational attraction g(E2, P) it exerts at the interior point P, which lies on E1, which now is just an aid to the calculations.
It turns out that g(E1, T(P)) and g(E2, P) have a very simple relationship, and this will allow us to use what we know about g for interior points to find its value for exterior points.
We can compute the x-component of g(E1, T(P)) by breaking up E1 into infinitesimal rectangular rods aligned with the x-axis, running the full length of the ellipsoid, with cross-sectional area dy dz. If we call the endpoints of the rod A and B, our formula for a rod gives us:
dgx(E1, T(P))
= G ρ dy dz [1/dist(T(P), A) – 1/dist(T(P), B)]
Similarly, we can compute the x-component of g(E2, P) by breaking up E2 into a matching set of rods with endpoints at T(A) and T(B) and a suitably rescaled cross-sectional area (αy αz)/(ay az) dy dz.
dgx(E2, P)
= G ρ (αy αz)/(ay az) dy dz [1/dist(P, T(A)) – 1/dist(P, T(B))]
= G ρ (αy αz)/(ay az) dy dz [1/dist(T(P), A) – 1/dist(T(P), B)]
= (αy αz)/(ay az) dgx(E1, T(P))
So for every such rod, the two quantities are related by the same factor, which is independent of the choice of rod, and when we integrate over all rods to get the totals they will be related in exactly the same way. What’s more, we get analogous relationships for the other coordinate directions.
gx(E2, P) = (αy αz)/(ay az) gx(E1, T(P))
gy(E2, P) = (αx αz)/(ax az) gy(E1, T(P))
gz(E2, P) = (αx αy)/(ax ay) gz(E1, T(P))
Since P is an interior point of the larger ellipsoid, E2, we already know the acceleration due to gravity here:
g(E2, P) = –4 π G ρ (Aα px, Bα py, Cα pz)
We have written these coefficients as Aα, Bα, Cα to emphasise that they need to be computed with the semi-axes of E2.
This lets us write the acceleration due to gravity from E1 at the exterior point T(P):
g(E1, T(P)) = –4 π G ρ [(ax ay az) / (αx αy αz)] (Aα T(P)x, Bα T(P)y, Cα T(P)z)
The effect of choosing a specific ellipsoid E1 here is less that it might seem. E1 singles out a family of confocal ellipsoids, but any other member of that family would do the same. Other than that, its semi-axes only appear as the product ax ay az, which is proportional to its volume, V = (4/3) π ax ay az, so we could remove this by specifying the mass of the ellipsoid, M, rather than its density and dimensions. The function T is defined with respect to E1 and E2, but T(P) could be any exterior point and we no longer need to relate it to P.
Acceleration due to gravity exterior to an ellipsoid E with uniform density and mass M
g(E, (x, y, z)) = [–3 M G / (αx αy αz)] (Aα x, Bα y, Cα z)
αx, αy and αz are the semi-axes of an ellipsoid confocal with E such that (x, y, z) lies on its surface.
The coefficients Aα, Bα, Cα are those described previously for the interior potential, computed with these α semi-axes.
So, for any of the members of a family of confocal ellipsoids with uniform density, the gravitational field, g, at a given exterior point will differ only in proportion to the mass of the ellipsoid.
As a consistency check, note that for a sphere A = B = C = ⅓, and αx = αy = αz = r, so the formula above becomes:
g(S, (x, y, z)) = [–M G / r3] (x, y, z)
which is just the inverse-square law, as expected. Also, if we choose a point on the boundary of the ellipsoid, the α semi-axes will be equal to the semi-axes of the ellipsoid itself, and the formula above will agree with the one for the interior.
What about the potential? Since the acceleration g at a given exterior point for different confocal ellipsoids only differs in proportion to their mass, the same must be true of the potential (given the constraint that it always goes to zero at infinity), so we can compute the potential by finding the value it would have if the ellipsoid was enlarged until its boundary contained the point in question, while its density was reduced in order to keep its mass unchanged.
Potential exterior to an ellipsoid E with uniform density and mass M
U(E, (x, y, z)) = [(3/2) M G / (αx αy αz)] (Aα x2 + Bα y2 + Cα z2 + Dα)
αx, αy and αz are the semi-axes of an ellipsoid confocal with E such that (x, y, z) lies on its surface.
The coefficients Aα, Bα, Cα, Dα are those described previously for the interior potential, computed with these α semi-axes.
The image on the left shows a slice through the surfaces of equal potential for a range of ellipsoids of revolution, all with the same density and volume. The axis of symmetry is vertical in this image, and the semi-axes vary from √2, √2, 1/2 (the most oblate) to √[2/3], √[2/3], 3/2 (the most prolate).
The curves of equal potential are perfect ellipses in the interior, but those in the exterior are not. It is only in the case of a sphere that the entire surface of the body has equal potential, and would be experienced as gravitationally level. For a prolate spheroid, where the polar axis is longer, the poles lie “uphill” from the equator, while for an oblate spheroid the equator is “uphill” from the poles. In both cases, the acceleration due to gravity, g, is weakest at the point that is farthest from the centre.
So far, we have considered only non-rotating bodies, which maintain their ellipsoidal shapes through the strength of the materials from which they are composed, which are assumed to be non-deformable solids. This is possible for small enough objects, but beyond a certain size the pressure due to gravity will be sufficient to make most materials flow, and so non-rotating bodies will take on a spherical shape.
What about rotating bodies? If a body undergoes uniform rotation with an angular velocity of ω around the z-axis, that adds a “centrifugal potential” of:
Ucen(x, y, z) = –½ω2(x2 + y2)
to the gravitational potential. The opposite of the gradient of this potential is:
acen = –∇Ucen = ω2 (x, y, 0)
|acen| = ω2 r
which is the centrifugal acceleration that material instantaneously at rest in the rotating frame will undergo, if subject to no other force.
If we combine this potential Ucen with the gravitational potential in the interior of an oblate spheroid, then the surfaces with a constant potential of Us will themselves be spheroids, satisfying the equation: