5  Newton’s shell theorem

import numpy as np
import sympy as sp
from sympy import pi, oo
from sympy.assumptions import assuming, Q
import matplotlib.pyplot as plt

We want to investigate spheres with an arbitrary radial density distribution \(\rho = \rho(r)\). The following results have been described 1687 in Newton’s Philosophiae Naturalis Principia Mathematica and have since then been known as Newton’s shell theorem.

Newton stated that:

Newton’s shell theorem

1. A spherical symmetric body causes a gravitational force on external objects as if all of its mass were concentrated in a point in its center.
2. If the body is a spherically symmetric shell, then no net gravitational force is exerted by the shell on any object inside regardless of the objects’s location within the shell.

Figure 5.1

The potential outside of a solid sphere of radius \(a\) and density law \(\rho=\rho(r)\) can be obtained by integration over spherical shells using the law of cosines: \[ V(r) = f \int_{r'=0}^a \int_{\theta=0}^\pi \int_{\phi=0}^{2\pi} \frac{\rho(r') {r'}^2 \, \mathrm d r' \sin\theta \, \mathrm d \theta \, \mathrm d \phi}{[{r'}^2 + r^2 -2 r r' \cos\theta]^{1/2}}. \]

Figure 5.2

We depart from the potential of one spherical shell of radius \(a\) with surface density \(\rho_F\), which can be evaluated as \[ V^F(r) = 2 \pi f \rho_F \int_{\theta=0}^\pi \frac{a^2 \sin\theta \, \mathrm d \theta}{\sqrt{r^2 + a^2 - 2ar\cos\theta}}. \]

The potential \(V^F(r)\) inside or outside of such a spherical shell with surface density \(\rho_F\) and radius \(a\) is \[ V^F(r) = 2 \pi f \frac{a}{r} \rho_F \left[ r + a - |r - a| \right]. \] Because \(r\) can be either less or greater than \(a\), we obtain two solution for \(V^F(r)\): \[ V^F(r) = \begin{cases} 4 \pi f \rho_F \frac{a^2}{r} & \quad\text{ if } r > a \\ 4 \pi f \rho_F a & \quad\text{ if } r < a \end{cases} \]

We now introduce a solid sphere with a density law \(\rho = \rho(r)\). To obtain the potential \(V(r)\) outside this sphere, we integrate

\[ V(r) = \int_{r'=0}^a V^F(r') \, \mathrm d r' = \frac{2 \pi f}{r } \int_{r'=0}^{a} \rho(r') [r + r' - |r - r'|] r' \, \mathrm d r' \]

To this end, we make use of SymPy, a Python library for symbolic mathematics.

We first define symbolic variables.

r, rp, a, f, rho, theta = sp.symbols("r r' a f rho theta", positive=True, real=True)

V = 2 * pi * f * rho * sp.integrate(a**2 * sp.sin(theta) / sp.sqrt(r**2 + a**2 - 2*a*r*sp.cos(theta)), (theta, 0, pi))
V.simplify()

\(\displaystyle \frac{2 \pi a f \rho \left(- \sqrt{a^{2} - 2 a r + r^{2}} + \sqrt{a^{2} + 2 a r + r^{2}}\right)}{r}\)

Since \[ \sqrt{r^2 + a^2 + 2 a r} = |r + a| = r + a \] and \[ \sqrt{r^2 + a^2 - 2 a r} = |r - a| = \begin{cases} r - a & \text{if } r > a \\ a - r & \text{if } r < a \end{cases} \] we have to be careful with the integrand.

Voutside = 2 * pi * f * rho / r * sp.integrate((r + rp - (r - rp)) * rp, (rp, 0, a))
Voutside

\(\displaystyle \frac{4 \pi a^{3} f \rho}{3 r}\)

To get the potential inside the sphere, i.e., for \(r < a\), we have to decompose the integral into two contributions.

First, we evaluate the potential due to a sphere of radius \(r\), to which we add the contribution of a hollow sphere with inner radius \(r\) and outer radius \(a\): \[ V(r) = \frac{4 \pi f}{r} \left[ \int_{r'=0}^r \rho(r') {r'}^2 \, \mathrm d r' + \int_r^a \rho(r') r r' \, \mathrm d r'\right] \]

In the interior of a homogeneous sphere with a mean density of \(\bar\rho\) we obtain for \(0 \le r \le a\) \[ V(r) = \frac{2 \pi f \bar\rho}{3} \left(3a^2 - r^2 \right) \]

Vinside = 4 * pi * f * rho / r * (sp.integrate(rp**2, (rp, 0, r)) + sp.integrate(rp * r, (rp, r, a)))
Vinside.simplify()

\(\displaystyle \frac{2 \pi f \rho \left(3 a^{2} - r^{2}\right)}{3}\)

5.1 Density laws

Now, we introduce non-trivial density laws, such as a linear density law \[ \rho(r) = \rho(0) - b r, \quad b > 0. \]

rho_m, rho_0, rho_R, R_E, b = sp.symbols("rho_m rho_0 rho_R R_E b", positive=True, real=True)

With the known mean density of the Earth we can evaluate Earth’s mass.

m_E = rho_m * 4 * pi / 3 * R_E**3
m_E

\(\displaystyle \frac{4 \pi R_{E}^{3} \rho_{m}}{3}\)

This is not very helpful. Instead, we can evaluate Earth’s mass using a volume integral over a density function, e.g., \[ m_E = 4 \pi \int_0^{R_E} (\rho_0 - br) r^2 \, \mathrm d r \] Using SymPy we obtain the expression

m_E_int = 4 * pi * sp.integrate((rho_0 - b * r) * r**2, (r, 0, R_E))
m_E_int.simplify()

\(\displaystyle \frac{\pi R_{E}^{3} \left(- 3 R_{E} b + 4 \rho_{0}\right)}{3}\)

This result still contains two unknown parameters \(b\) and \(\rho_0\).

We exploit the fact that there must hold \[ \rho(R_E) = \rho(0) - b R_E. \]

The density at the Earth’s surface is about \(2700\) \(kg/m^3\).

eq1 = m_E_int - m_E
eq2 = rho_0 - b * R_E - rho_R

coeffs = sp.solve((eq1, eq2), (rho_0, b),dict=True)
coeff_rho_0 = coeffs[0][rho_0]
coeff_b = coeffs[0][b]

coeff_rho_0

\(\displaystyle - 3 \rho_{R} + 4 \rho_{m}\)

coeff_b

\(\displaystyle \frac{- 4 \rho_{R} + 4 \rho_{m}}{R_{E}}\)

rho_num = coeff_rho_0.subs(rho_R, 2700).subs(rho_m, 5515)

b_num = coeff_b.subs(rho_R, 2700).subs(rho_m, 5515).subs(R_E, 6.371e6)

5.2 Gravitational attraction inside the Earth

It is known that the gravitational potential inside a sphere with a radially symmetric density function is \[ V(r) = \frac{4 \pi f}{r} \left[ \int_{r'=0}^r \rho(r') {r'}^2 \, \mathrm d r' + \int_r^{R_E} \rho(r') r r' \, \mathrm d r'\right] \] The first integral amounts for the potential at the surface of a sphere with radius , whereas the second computes the potential taken at the inner surface of a hollow sphere.

With the constants derived above we want to evaluate these integrals.

i1 = 4 * pi * f / r * sp.integrate((rho_0 - b * rp) * rp**2, (rp, 0, r))
i1.simplify()

\(\displaystyle \frac{\pi f r^{2} \left(- 3 b r + 4 \rho_{0}\right)}{3}\)

i2 = 4 * pi * f * sp.integrate((rho_0 - b * rp) * rp, (rp, r, R_E))
i2.simplify()

\(\displaystyle \frac{2 \pi f \left(- 2 R_{E}^{3} b + 3 R_{E}^{2} \rho_{0} + 2 b r^{3} - 3 r^{2} \rho_{0}\right)}{3}\)

Next, we add both contributions to obtain the potential inside the Earth.

#V = sp.Function("V")
V = i1 + i2

V.simplify()

\(\displaystyle \frac{\pi f \left(- 4 R_{E}^{3} b + 6 R_{E}^{2} \rho_{0} + b r^{3} - 2 r^{2} \rho_{0}\right)}{3}\)

We obtain the gravitational attraction by taking the gradient of \(V\) in the direction of \(\mathbf r\):

\[ g(r) = \nabla V(r) \cdot \hat{\mathbf r} \]

g = sp.diff(V, r)
g.simplify()

\(\displaystyle \frac{\pi f r \left(3 b r - 4 \rho_{0}\right)}{3}\)

5.3 Visualization of the result

We plot the curve of \(g(r)\) inside the Earth.

Further, we check if \(g(r)\) has an extremum inside the Earth.

f_ = 6.674e-11
g_num = lambda r: np.pi * f_ * r * (3 * b_num * r - 4 * rho_num) / 3

rstar = sp.solve(sp.diff(g, r), (r))
rstar[0]

\(\displaystyle \frac{2 \rho_{0}}{3 b}\)

rextrem = rstar[0].subs(rho_0, rho_num).subs(b, b_num)
rextrem / 1e6

\(\displaystyle 5.26578804026051\)

6.371e6 - rextrem

\(\displaystyle 1105211.95973949\)

The extremal value of \(g(r)\) in the interior of the Earth appears at \(r=5.26 \times 10^6\) m, i.e., at a depth of 1105 km.

rr = np.linspace(0, 6.367e6, 41)

plt.plot(rr, g_num(rr))
plt.scatter(rextrem, g_num(rextrem), color='r')
plt.grid(True)
plt.xlabel("r in m")
plt.ylabel("g(r)");

5.4 Moment of inertia

There are two constraints for density models: The mass of the Earth, and its moment of inertia. Both are moments of the density distribution. The mass is the second moment of the radial density distribution, whereas the mean moment of inertia is the scaled fourth moment of the radial density distribution. The moment of inertia can be inferred from a thin spherical shell of radius \(a\),

\[ I = \frac{2}{3} m a^2. \]

The contribution of a small mass \(\mathrm d m\) is therefore \[ \mathrm d I = \frac{2}{3} a^2 \, \mathrm d m. \] The total moment of intertia an be calculated with the integral \[ I=\frac{2}{3} \int\limits_0^{R_E} \rho(r) r^2 \, \mathrm d V = \frac{8\pi}{3} \int\limits_0^{R_E} \rho(r) r^4 \, \mathrm d r, \] with the help of which we compute the moment of inertia factor \(\alpha\), a dimensionless number, where \[ \alpha = \frac{I}{mR_E^2}. \] This factor is \(0.4\) for a solid sphere.

Check using SymPy:

(8 * pi / 3 * 
    sp.integrate(rho_m * r**4, (r, 0, R_E)) / 
    m_E / R_E**2
    ).subs(rho_m, 5515).subs(R_E, 6.371e6)

\(\displaystyle \frac{2}{5}\)

\(\alpha\) should be smaller for a sphere when there is an increasing density with depth.

For our density model, we observe \(\alpha=0.3319\), which is slightly more than reported in the literature, where \(\alpha=0.3307144\) (Williams, 1994).

Check:

(8 * pi / 3 * 
    sp.integrate((rho_num - b_num * r) * r**4, (r, 0, R_E)) / 
    m_E / R_E**2
    ).subs(rho_m, 5515).subs(R_E, 6.371e6)

\(\displaystyle 0.331943185252342\)