2.2. Spherical systems: Newton’s shell theorems¶
For spherical mass distributions, Newton proved two fundamental theorems that significantly simplify all work with spherical mass distributions and, in particular, that of solving the Poisson equation. These are:
Newton’s first shell theorem: A body that is inside a uniform spherical shell of matter experiences no net gravitational force from that shell.
Newton’s second shell theorem: The gravitational force on a body that lies outside a uniform spherical shell of matter is the same as it would be if all of the shell’s matter were concentrated into a point at its center.
A direct, mathematical proof is most easily based on Gauss’s law, which is a direct consequence of the Poisson equation. Let’s integrate the Poisson equation over an arbitrary volume \(V\) containing a total mass \(M\) and bounded by the surface \(S\) and use the divergence theorem (Equation B.5) to turn the volume integral into a surface integral. Then we obtain that \begin{align} 4\pi\,G\,M & = 4\pi\,G\,\int_V\mathrm{d}V\,\rho = \int_V\mathrm{d}V\,\nabla^2 \Phi\\ & = \int_S\mathrm{d}S\,\left(\vec{\hat{n}}\cdot\nabla \Phi\right) = -\int_S\mathrm{d}S\,\left(\vec{\hat{n}}\cdot\vec{g}\right)\,, \end{align} where \(\vec{\hat{n}}\) is the unit vector perpendicular to the surface \(S\). Thus, the integral of the component of the gravitational field perpendicular to a given surface is equal to the mass contained within that surface (multiplied by \(-4\pi G\)).
To prove Newton’s first shell theorem, we consider a spherical shell \(S_a\) centered on the origin with radius \(a\) and we integrate the gravitational field over a similar spherical shell \(S_b\) with radius \(b < a\). By symmetry, the gravitational field can only depend on \(r\) and for a non-zero field would therefore be constant \(g_b\) on \(S_b\); its integral over \(S_b\) is therefore \(4\pi b^2 \,g_b\)—the surface of the shell times the gravitational field. By Gauss’ law, this should equal \(-4\pi\,G\) times the enclosed mass, which is zero because \(S_b\) is within \(S_a\). Therefore, the gravitational field \(g_b = 0\). This holds for all \(b < a\), which proves Newton’s first shell theorem.
Newton’s second shell theorem can be proven in a similar way. We now integrate over a shell \(S_c\) with \(c > a\). Again, the gravitational field on this shell is a constant \(g_c\) and the integral is equal to \(4\pi c^2 g_c\). Because \(S_c\) is outside of \(S_a\), all of \(S_a\)’s mass is contained within \(S_c\) and Gauss’s law now implies that this should equal \(-4\pi\,G\) times the mass of the shell \(S_a\). Because this equality does not depend on the radius \(a\) of \(S_a\), it would be the same if we shrunk the shell to a point, thus proving Newton’s second shell theorem.
Newton originally proved his first shell theorem using more geometric means that are more intuitive. We consider the force on a point interior to a thin shell from a small section of the shell obtained from the two intersections of a cone centered on this point with a small opening angle as in Figure 2.2.
[8]:
figure(figsize=(6,6))
gca().add_patch(Circle((0.2,0.1),radius=1,fc='none',ec='k',zorder=0))
gca().add_patch(Circle((0.2,0.1),radius=1.05,fc='none',ec='k',zorder=0))
gca().add_patch(Circle((0.5,0.6),radius=.04,fc='k',ec='none',zorder=2))
plot((-1.,1.1,),(5/7*-1-5/14+6/10,5/7*1.1-5/14+6/10),zorder=1,color='#1f77b4')
plot((-.6,0.9,),(6/4*-.6-6/8+6/10,1.2),zorder=1,color='#1f77b4')
# Find intersection points
m, b= 5/7.,-5/14+6/10
pp= numpy.roots((1+m**2.,-(0.4-2*m*(b-0.1)),0.04-1.025**2+(b-0.1)**2.))
p= (pp[pp<0.2],m*pp[pp<0.2]+b)
q= (pp[pp>0.2],m*pp[pp>0.2]+b)
m, b= 6/4.,-6/8+6/10
pp= numpy.roots((1+m**2.,-(0.4-2*m*(b-0.1)),0.04-1.025**2+(b-0.1)**2.))
qp= (pp[pp>0.2],m*pp[pp>0.2]+b)
pp= (pp[pp<0.2],m*pp[pp<0.2]+b)
# intersection lines
m,b= (pp[1]-p[1])/(pp[0]-p[0]), -(pp[1]-p[1])/(pp[0]-p[0])*p[0]+p[1]
#m,b= -(p[0]-0.2)/(p[1]-0.1),p[1]+(p[0]-0.2)/(p[1]-0.1)*p[0]
plot((p[0]+0.525,p[0]-0.2),(m*(p[0]+0.525)+b,m*(p[0]-0.2)+b),color='#ff7f0e')
m,b= (qp[1]-q[1])/(qp[0]-q[0]), -(qp[1]-q[1])/(qp[0]-q[0])*q[0]+q[1]
#m,b= -(q[0]-0.2)/(q[1]-0.1),q[1]+(q[0]-0.2)/(q[1]-0.1)*q[0]
plot((q[0]+0.2,q[0]-0.3),(m*(q[0]+0.2)+b,m*(q[0]-0.3)+b),color='#ff7f0e')
plot((p[0],pp[0],q[0],qp[0]),(p[1],pp[1],q[1],qp[1]),'ko',ms=5.)
text(-0.5,0.,r'$d_1$',size=18.)
text(0.5,0.825,r'$d_2$',size=18.)
text(-1.225,-0.8,r'$\delta m_1 \propto d_1^2$',size=18.)
text(0.9,1.1,r'$\delta m_2 \propto d_2^2$',size=18.)
xlim(-1.25,1.5)
ylim(-1.25,1.5)
gca()._frameon= False
gca().xaxis.set_visible(False)
gca().yaxis.set_visible(False)
Figure 2.2: Newton’s first shell theorem.
For a narrow cone, the surface of the intersection between the cone and the shell is approximately that of the cone and a plane; such an intersection is a conical section and in this case in particular, an ellipse. The major and minor axis are the segment of the orange line between the two dots in Figure 2.2 (major or minor depending on the projection). For small opening angles, the ratio of the lengths of these axes is equal to the ratio of the distances between the center of the cone and the intersections: \(d_1/d_2\) in the figure. Therefore, the ratio of the areas of the ellipses is \(d_1^2/d_2^2\) and, because the shell has uniform density and thickness, the ratio of the masses is the same. Because of the \(1/r^2\) dependence of the gravitational force, the gravitational forces from the two intersections are then equal in magnitude and opposite in sign and they therefore cancel. This holds for any narrow cone centered on any point within the shell and, therefore, there is no net force from the shell on any interior point.
The first theorem implies that for a spherical mass distribution, mass outside of the current radius of a body has no influence on the motion of that body (at the present time). In particular, for a body on a circular orbit, the mass outside of this circle has no effect on the entire orbit. We will see in Chapter 7 that this is not the case for flattened mass distributions (the geometric proof of Newton’s first theorem immediately makes it clear why it does not hold for flattened distributions).
Newton’s second shell theorem implies that the gravitational potential outside of a shell of radius \(R\) is \begin{equation}\label{eq-spherpot-outside} \Phi_{\mathrm{shell}}(r>R) = -\frac{GM_{\mathrm{shell}}}{r}\,. \end{equation} The first theorem together with the requirement that the potential is continuous at \(R\) then implies that the gravitational potential within a shell of radius \(R\) is equal to \begin{equation} \Phi_{\mathrm{shell}}(r < R) = -\frac{GM_{\mathrm{shell}}}{R}\,. \end{equation} An example of the potential of a shell at \(r=1.5\) is displayed in Figure 2.3.
[6]:
from galpy.potential import SphericalShellPotential
sp= SphericalShellPotential(a=1.5)
rs= numpy.linspace(0.1,10.,101)
galpy_plot.plot(rs,[sp(r,0.) for r in rs],
yrange=[-.8,0.],
xlabel=r'$r$',ylabel=r'$\Phi_{\mathrm{shell}}$');
Figure 2.3: The potential of a spherical shell.
To obtain the gravitational potential at radius \(r\) due to a spherical mass distribution \(\rho(r')\), we can therefore sum the contributions from all shells with mass \(\mathrm{d}M(r') = 4\pi G\rho(r')r'^2\mathrm{d}r'\) inside and outside of \(r\) as follows \begin{equation}\label{eq-spherpot} \Phi(r) = -4\pi\,G\,\left[\frac{1}{r}\int_0^r\mathrm{d}r'\,\rho(r')\,r'^2+\int_r^\infty\mathrm{d}r'\,\rho(r')\,r'\right]\,. \end{equation} One can show that this satisfies the Poisson equation using the Laplacian in spherical coordinates (Equation A.8). Thus, for spherical mass distributions, we can compute the potential \(\Phi\) using a simple quadrature for any mass distribution \(\rho(r)\).
As an example, we discretize the mass distribution of a Plummer profile (see below) into a small number of shells where each shell contains the mass between the radius of the shell just inside it and the shell’s radius (and the innermost shell contains all mass up to its radius). We can do this with the following function, which creates a galpy potential that is a set of discrete shells defined in this way for any spherical potential:
[7]:
from scipy import integrate
from galpy.potential import SphericalShellPotential
def discretize_potential_into_shells(pot,rmin,rmax,dr):
rs= numpy.arange(rmin,rmax+dr,dr)
dM= numpy.empty_like(rs)
dM[0]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
0.,rs[0])[0]
for ii in range(1,len(rs)):
dM[ii]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
rs[ii-1],rs[ii])[0]
return [SphericalShellPotential(amp=dm,a=r) for (dm,r) in zip(dM,rs)]
Discretizing a Plummer profile in this way with scale parameter \(b=1.5\), using shells spaced \(\Delta r = 1\) apart out to \(r = 3\times b\) gives the potential shown in the left panel of Figure 2.4.
[8]:
from galpy.potential import PlummerPotential, evaluatePotentials
pp= PlummerPotential(amp=1.,b=1.5)
rs= numpy.linspace(0.1,10.,101)
figure(figsize=(12,3.45))
subplot(1,3,1)
rmin, rmax, dr= 0.5,4.5,1.
discrete_pp= discretize_potential_into_shells(pp,rmin,rmax,dr)
line_approx= galpy_plot.plot(rs,[evaluatePotentials(discrete_pp,r,0.)
for r in rs],
yrange=[-.75,0.],gcf=True,
xlabel=r'$r$',ylabel=r'$\Phi$',
label=r'$\mathrm{Shell\ approximation}$'\
+'\n'+r'$r_\mathrm{max}=3b, \Delta r=1$')
line_plummer= galpy_plot.plot(rs,pp(rs,0.),overplot=True,zorder=0,
label=r'$\mathrm{Plummer}$')
legend(handles=[line_approx[0],line_plummer[0]],
fontsize=17.,loc='lower right',frameon=False)
subplot(1,3,2)
rmin, rmax, dr= 0.5,4.5,0.1
discrete_pp= discretize_potential_into_shells(pp,rmin,rmax,dr)
line_approx= galpy_plot.plot(rs,[evaluatePotentials(discrete_pp,r,0.)
for r in rs],
yrange=[-.75,0.],gcf=True,
xlabel=r'$r$',
label=r'$\mathrm{Shell\ approximation}$'\
+'\n'+r'$r_\mathrm{max}=3b, \Delta r=0.1$')
line_plummer= galpy_plot.plot(rs,pp(rs,0.),overplot=True,zorder=0,
label=r'$\mathrm{Plummer}$')
legend(handles=[line_approx[0],line_plummer[0]],
fontsize=17.,loc='lower right',frameon=False)
gca().tick_params(labelleft=False)
subplot(1,3,3)
def discretize_potential_into_shells_addinf(pot,rmin,rmax,dr):
rs= numpy.arange(rmin,rmax+dr,dr)
dM= numpy.empty_like(rs)
dM[0]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
0.,rs[0])[0]
for ii in range(1,len(rs)-1):
dM[ii]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
rs[ii-1],rs[ii])[0]
dM[-1]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
rs[ii],numpy.inf)[0]
return [SphericalShellPotential(amp=dm,a=r) for (dm,r) in zip(dM,rs)]
rmin, rmax, dr= 0.5,7.5,0.1
discrete_pp= discretize_potential_into_shells_addinf(pp,rmin,rmax,dr)
line_approx= galpy_plot.plot(rs,[evaluatePotentials(discrete_pp,r,0.)
for r in rs],
yrange=[-.75,0.],gcf=True,
xlabel=r'$r$',
label=r'$\mathrm{Shell\ approximation}$'\
+'\n'+r'$r_\mathrm{max}=5b, \Delta r=0.1$')
line_plummer= galpy_plot.plot(rs,pp(rs,0.),overplot=True,zorder=0,
label=r'$\mathrm{Plummer}$')
legend(handles=[line_approx[0],line_plummer[0]],
fontsize=17.,loc='lower right',frameon=False)
gca().tick_params(labelleft=False)
tight_layout(w_pad=-1.1);
Figure 2.4: Shell approximations of a Plummer potential.
The orange curve shows the actual potential for a Plummer profile (given in Section 2.4.3 below). We see the discreteness of the approximation clearly at \(r \lesssim 2b\), but at larger \(r\) the approximation becomes largely smooth because most of the mass is within a few times \(b\) for a Plummer profile. We also notice that the approximation, while smooth, remains too high at large \(r\). Decreasing the spacing of the shells improves the agreement at small \(r\) as shown in the middle panel of Figure 2.4. But there is still an overall offset between the shells and the true potential. This is because we are ignoring the mass at \(r > 3b\), thus clearly showing that the potential—unlike the force—at \(r\) depends on the mass profile outside of \(r\). We therefore edit our shell-discretization function to let the final shell contain all mass out to infinity:
[13]:
def discretize_potential_into_shells_addinf(pot,rmin,rmax,dr):
rs= numpy.arange(rmin,rmax+dr,dr)
dM= numpy.empty_like(rs)
dM[0]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
0.,rs[0])[0]
for ii in range(1,len(rs)-1):
dM[ii]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
rs[ii-1],rs[ii])[0]
dM[-1]= 4.*numpy.pi*integrate.quad(lambda r: r**2*pot.dens(r,0.),
rs[ii],numpy.inf)[0]
return [SphericalShellPotential(amp=dm,a=r) for (dm,r) in zip(dM,rs)]
Using this function with shells out to \(r = 5b\) then gives the result from the right panel in Figure 2.4. Now we have almost perfect agreement for all but the innermost part of the potential well.
The only non-zero component of the gravitational field is the radial component \(g_r\), which is given by \begin{equation}\label{eq-sphere-gr} g_r(r) = -\frac{G\,M(< r)}{r^2}\,, \end{equation} where \(M(<r)\) is the mass contained within radius \(r\). This can be derived directly from Newton’s theorems or by taking the derivative from Equation (2.21) for the potential. Because \(g_r(r) = -\mathrm{d} \Phi(r) / \mathrm{d} r\), this gives an alternate relation between the enclosed mass and the potential \begin{equation}\label{eq-spherpot-encmass} \Phi(r) = -G\,\int_r^\infty \mathrm{d}r'\,\frac{M(< r')}{r'^2} = -4\pi\,G\, \int_r^\infty \mathrm{d}r'\,\frac{1}{r'^2}\,\int_0^{r'}\mathrm{d}r''\,r''^2\,\rho(r'') \,. \end{equation} You can show that Equation (2.23) is equal to Equation (2.21) by changing the order of the integration (carefully). These forms can make the calculation of the potential easier if the enclosed mass has a known, simple form.
.
.
.