3. Gravitation¶

In this chapter, we introduce the basics of gravitational dynamics: how does mass give rise to forces that lead to the motion of planets, stars, gas, and dark matter in the Universe? The next chapters discuss how masses move under the influence of the gravitational force—these motions we refer to as “orbits”—which is the realm of classical mechanics, and we introduce the important concept of dynamical equilibrium.

In this first part, we apply these concepts to mass distributions that are spherically symmetric. Real galaxies look like this:

(Credit: M101: European Space Agency & NASA; NGC 660: Gemini Observatory, AURA)

One might therefore wonder how useful it is to study spherical mass distributions, when most galaxies are so far from being spherically symmetric. We will see in later chapters that, even though mass distributions can be quite far from spherical, many of the properties of non-spherical distributions can be understood by approximating these distributions with some equivalent spherical distribution. Similarly, many of the concepts that we will introduce in this chapter (like that of dynamical time and the circular velocity) remain similar for non-spherical distributions. And spherical mass distributions provide a simple manner to introduce many of the concepts that are relevant for all galaxies.

3.1. Matter and the gravitational field¶

A theory for the motion of objects under the influence of gravity requires two ingredients: how mass gives rise to a gravitational force and how objects move under the influence of this force. In classical physics, the laws for each of these go back to Newton and contemporaries: Newton’s law of universal gravitation that gives the force \(F\) from a mass \(M\) acting on a second mass \(m\) in terms of the distance \(r\) between the two as \(F = -GMm/r^2\), and Newton’s second law that gives the relation between force and acceleration \(a\) for an object with mass \(m\) as \(F = ma\). A more modern understanding of the law of universal gravitation, however, emphasizes that the three-dimensional force \(\vec{F}(\vec{x})\) derives from a scalar-valued gravitational potential \(\Phi(\vec{x})\) through:

\begin{equation}\label{eq-force-gradient-potential} \vec{F}(\vec{x}) = -m\,\nabla \Phi(\vec{x})\,, \end{equation}

and replaces Newton’s inverse-square law of gravitation with the Poisson equation as the fundamental equation relating mass and gravitational force:

\begin{equation}\label{eq-sphergrav-poisson} \nabla^2 \Phi(\vec{x}) = 4\pi\,G\,\rho(\vec{x})\,, \end{equation}

where \(\rho(\vec{x})\) is the mass density, \(G\) is the gravitational constant, and \(\nabla^2\) is the Laplace operator.

The reason that the Poisson equation is more properly considered to be the fundamental equation between mass and gravitational force is that it is the direct Newtonian limit of Einstein’s field equation \(R_{\mu\nu} -R\,g_{\mu\nu}/2 = 8\pi G\,T_{\mu\nu}\) in the general theory of relativity, which is our best theory of gravity so far. This is not a book about the general theory of relativity and the origin of Einstein’s field equation is typically not important for understanding galaxy dynamics, but a key point is that Einstein’s field equation reduces to the Poisson equation in the limit that velocities \(v\) and the gravitational potential are small compared to the speed of light \(c\) (the potential \(\Phi\) has units of velocity squared, so this limit is \(|\Phi|/c^2 \ll 1\) and \(v/c\ll 1\)). This limit always applies on the scale of galaxies, with the notable exception of gravitational lensing by galaxies as discussed in Chapter 16. For the interested reader, Appendix C has a self-contained discussion of Einstein’s field equation and how it reduces to the Poisson equation in the low-velocity limit.

Before continuing our discussion, it is worth noting that the Poisson equation (or its generalization, Einstein’s field equation) is a hypothesis. Like any other physical theory, we test the Poisson equation by making predictions derived using this equation and testing these predictions against observational data. Gravity as defined by the Poisson equation is extraordinarily well tested using laboratory experiments and solar-system gravitational dynamics: it is known to hold to one-part-in-one-billion in these settings (Adelberger et al. 2003; Will 2014). Partly because of these tests, most astrophysicists assume that the Poisson equation holds on the scales of galaxies and on the scale of the Universe (as Einstein’s equation) and that any anomalies that may result from its application are due to new forms of mass or energy that enter the mass density \(\rho\) or the energy-momentum tensor \(T_{\mu\nu}\), rather than indicating a problem with the Poisson equation. Large anomalies are indeed known to exist: these are dark matter, which is a \(\approx100\%\) anomaly on the scales of galaxies, and dark energy, which is a \(\approx100\%\) anomaly on the largest cosmological scales (meaning that ignoring them leads to \(\mathcal{O}(1)\) discrepancies with the observations).

Throughout this book, we will largely assume that the Poisson equation holds and, thus, that dark matter (and dark energy, but it is less important for our purposes) is a new form of matter whose density distribution can be studied using the Poisson equation. Theories that attempt to account for the anomalies typically interpreted as dark matter and dark energy without introducing new forms of matter and energy do so by modifying the Poisson equation or Einstein’s equation. Such modifications have to account for the fact that the Poisson equation holds to high accuracy in the solar system and, thus, can only change the Poisson equation in vastly different physical regimes.

Starting from the Poisson equation, we can derive Newton’s inverse-square law by computing the gravitational potential for a point mass \(M\) at position \(\vec{x}\). Because the Laplacian is a differential operator, we can move the origin of the coordinate system to any position, and calculations are most straightforward if we position the mass at the origin. The density in this case is \(\propto\delta(\vec{x})\), where \(\delta(\cdot)\) is Dirac’s delta function (Chapter B.1.1) . To show that the potential \(\Phi(\vec{x}) = -GM/r\), where \(r = |\vec{x}|\), is the solution of the Poisson equation in this case, we compute the Laplacian of \(-GM/r\) for \(r \neq 0\) using the expression for the Laplacian in spherical coordinates (Equation A.12):

\begin{align} \nabla^2 \left(\frac{-GM}{r}\right) & = -GM\,\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\,\frac{\partial}{\partial r}\left[\frac{1}{r}\right]\right)\\ & = GM\,\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\,\frac{1}{r^2}\right)\\ & = 0\,. \end{align}

Thus, the density is zero for all \(r\neq 0\). Finally, to show that at \(r=0\), the density corresponding to \(\Phi(\vec{x}) = -GM/r\) is the correct delta function, we integrate the Laplacian of \(\Phi/(4\pi G)\) over a small spherical volume \(V\) of radius \(R\) with surface \(S\), because this should equal the mass \(M\):

\begin{align} \int_V \mathrm{d}V\, \nabla^2\left(-\frac{GM}{4\pi G\,r}\right) & = -\frac{M}{4\pi}\,\int_V \mathrm{d} V\, \nabla\cdot\left(\nabla\frac{1}{r}\right)\\ & = -\frac{M}{4\pi}\,\int_S \mathrm{d} S\, \vec{\hat{r}}\cdot \nabla\frac{1}{r}\\ & = -\frac{M}{4\pi}\,\int_S \mathrm{d} S\, \vec{\hat{r}}\cdot \vec{\hat{r}}\,\frac{\partial}{\partial r}\left(\frac{1}{r}\right)\\ & = \frac{M}{4\pi}\,\int_S \mathrm{d} S\, \frac{1}{r^2}\\ & = 4\pi\,\frac{M}{4\pi}\,\frac{R^2}{R^2}\\ & = M\,. \end{align}

Here, we have used the divergence theorem in going to the second line and the expression for the gradient in spherical coordinates in going to the third line (Equation A.11); \(\vec{\hat{r}}\) is the unit vector in the radial direction, which is perpendicular to the surface \(S\). Thus, the gravitational potential of a point mass \(M\) at distance \(r\) is \(\Phi = -GM/r\). Using that the gravitational force is the gradient of the potential from Equation \(\eqref{eq-force-gradient-potential}\), we find Newton’s law of gravity: at a three-dimensional position \(\vec{x}\) a distance \(r\) from the point-mass \(m\)

\begin{equation} \vec{F}(\vec{x}) = -\frac{GMm}{r^2}\,\vec{\hat{x}}\,. \end{equation}

Putting the position \(\vec{x}_0\) of the mass \(M\) back in explicitly, this force law becomes

\begin{equation} \vec{F}(\vec{x}) = -\frac{GMm}{|\vec{x}-\vec{x}_0|^3}\,(\vec{x}-\vec{x}_0)\,, \end{equation}

where we have used the standard simplification that the unit vector along \((\vec{x}-\vec{x}_0)\) is \((\vec{x}-\vec{x}_0)/|\vec{x}-\vec{x}_0|\). Because the force falls off as \(1/r^2\) rather than following something like a fast exponential decline, it is called a long-range force. In particular, if a point mass \(m\) is surrounded by a density of point masses that is uniform, then the amount of mass \(M\) in a shell at distance \(R\) is \(\propto R^2\), which combined with the \(1/R^2\) behavior of the force means that the force acting on \(m\) from shells at different \(R\) is approximately constant. In a constant-density medium, there are many more such shells at large distances than at small distances, and therefore the force is dominated by the total contribution from distant shells, rather than that from a few nearby shells.

A pair of functions (\(\Phi,\rho\)) that solve the Poisson equation \((\nabla^2 \Phi = 4\pi G \rho)\) is known as a potential-density pair. Because the Laplacian is a linear differential operator, we have that a linear combination of solutions to the Poisson equation is itself a solution: for potential-density pairs (\(\Phi_1,\rho_1\)) and (\(\Phi_2,\rho_2\)), the pair (\(\Phi_1+\Phi_2,\rho_1+\rho_2\)) is also a solution. A consequence of this is that the gravitational potential for a set of \(N\) point masses is simply given by the sum of the potentials for the individual point masses: If an object at \(\vec{x}\) with mass \(m\) is at a distance \(d_i = |\vec{x}-\vec{x}_i|\) from point masses \(M_i\) at positions \(\vec{x}_i\), the total gravitational potential is

\begin{equation} \Phi(\vec{x}) = \sum_i \frac{-GM_i}{d_i}\,. \end{equation}

Similarly, the gravitational force is

\begin{equation} \vec{F}(\vec{x}) = \sum_i \frac{-GM_i m}{d_i^3}\,(\vec{x}-\vec{x_i})\,. \end{equation}

Note that this follows directly from the Poisson equation, while deriving it from Newton’s law of universal gravitation would require the additional assumption that forces add up linearly (which is of course baked into the Poisson equation).

The mass of galaxies is contained in discrete chunks, whether they be stars, putative dark-matter particles, or the atoms and molecules of the interstellar medium. Even though this matter is discrete, the overall distribution of mass is rather uniform, and the gravitational force even between large chunks like stars is therefore dominated by distant bodies (see the argument above). Therefore, we can approximate the density in a galaxy as a smooth function, rather than as a sum over discrete bodies. From the Poisson equation, this means that the gravitational potential and gravitational force are smooth functions as well: Because the density is a second derivative of the potential, the potential is essentially a double integral of the density and, therefore, much smoother than the density.

Newton’s second law (which will be discussed in more detail in the next chapter), states that mass times acceleration equals force. The mass that appears in this equation is the same mass as that appears in the equation between the gravitational potential and force (Equation [\(\ref{eq-force-gradient-potential}\)]; or in Newton’s law of gravitation if you wish). Therefore, we have that

\begin{equation} F = -m\,\nabla \Phi(\vec{x}) = m\,\vec{a}\,, \end{equation}

or

\begin{equation} -\nabla \Phi(\vec{x}) = \vec{a}\,. \end{equation}

The motion of an object in a smooth, external gravitational potential therefore does not depend on its mass. This fails when the field is not external, that is, when the object’s mass has an effect on its surrounding mass distribution, which in turn affects its motion through Newton’s second law. But for many applications of galaxy dynamics, a smooth, external gravitational potential is an excellent approximation. It therefore makes sense to introduce the gravitational field \(\vec{g}(\vec{x})\)—the force per unit mass— as

\begin{equation} \vec{F}(\vec{x}) = m\vec{g}(\vec{x})\,. \end{equation}

because then in a smooth, external potential we have that

\begin{equation}\label{eq-accel-is-field} \vec{a} = \vec{g}\,. \end{equation}

This is known as the weak equivalence principle or the universality of free fall: all objects fall the same in an external gravitational field, whether they be feathers, stones, or stars. Einstein made this principle the centerpiece of his theory of relativity.

Because most of the time we do not need to consider an object’s mass to discuss its motion under gravity, we typically deal only with the gravitational field and much of the literature on galactic dynamics uses the terms “gravitational force” and “gravitational field” interchangeably and often uses the force symbol \(\vec{F}\) when really the field \(\vec{g}\) is meant. In this book, I will attempt to correctly use the terms “force” and “field”. Similarly, because of the weak equivalence principle, the gravitational field and the acceleration caused by it are the same, and the terms acceleration and force/field are often used interchangeably. Obviously we can only do this when the force is that due to gravity.

When working with physical quantities and equations, it is often useful to keep their units in mind. From the universality of free fall, we know that the units of gravitational field are the same as those of acceleration: length over time squared. The gravitational field is the spatial derivative of the gravitational potential. Therefore, gravitational potential has units of length squared over time squared, or more simply, velocity squared. Conversely, from the Poisson equation, we have that \(G\times\mathrm{density}\) has the same units as the spatial derivative of the gravitational field: units of inverse-time squared. \(G\) itself therefore has units of length cubed over mass over time squared. For galactic systems, \(G\) is most usefully expressed as (using the CODATA 2018 version of the recommended values of the fundamental physical constants)

\begin{equation} G = 4.301 \times 10^{-3}\,\mathrm{pc}\,\left(\mathrm{km/s}\right)^{2}\,M_\odot^{-1}\,. \end{equation}

where \(M_\odot\) is the mass of the Sun. We can also write this as

\begin{equation} G = 4.301\,\mathrm{kpc}\,\left(100\,\mathrm{km/s}\right)^{2}\,\left(10^{10}\,M_\odot\right)^{-1}\,. \end{equation}

Because gravity is such a weak force, measuring \(G\) precisely is difficult and the relative uncertainty of the current measurement is \(2\times 10^{-5}\). This is a much larger relative uncertainty than those of other physical constants, which have typical relative uncertainties of \(\approx 10^{-10}\). However, essentially all of astrophysics is only sensitive to the combination \(GM_\odot\), which is known to a relative uncertainty of \(\approx 10^{-10}\), because it is measured using the orbits of spacecraft in the solar system. It is worth remembering that whenever you see a quoted measured mass outside of the solar system, what is really determined is \(G\) times that mass.

To sum up: The Poisson equation is the fundamental equation one has to solve to obtain the gravitational force due to any mass distribution. Because of this fundamental relation, we will use the terms mass distribution and (gravitational) potential interchangeably (where “gravitational” is typically implied in this context if it is not mentioned in front of “potential”). The (negative) gradient of the potential gives the gravitational field that gives rise to motion and is therefore the only quantity that has physical significance. Thus, we can add or subtract any constant from the potential without changing the dynamics; whenever possible we shall fix this constant such that the potential equals zero at \(r = \infty\).

3.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 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 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 theorem, which is a direct consequence of the Poisson equation. Let’s integrate the Poisson equation over an arbitrary volume \(V\) containing total mass \(M\) and bounded by the surface \(S\) and use the divergence theorem 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\\ & = \phantom{4\pi\,G\,}\int_V\mathrm{d}V\,\nabla^2 \Phi\\ & = \phantom{4\pi\,G\,}\int_S\mathrm{d}S\,\left(\vec{\hat{n}}\cdot\nabla \Phi\right)\\ & = \phantom{\,G\,}-\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’ theorem, 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’ theorem 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 the following figure:

[6]:

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)

../_images/chapters_I-01.-Potential-Theory-and-Spherical-Mass-Distributions_10_0.png

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 the figure above (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 above. 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 8 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:

[8]:

from galpy.potential import SphericalShellPotential
from galpy.util import plot as galpy_plot
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}}(r)$');

../_images/chapters_I-01.-Potential-Theory-and-Spherical-Mass-Distributions_14_0.png

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.12). 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:

[9]:

from scipy import integrate
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)]

Applying this for a Plummer profile with scale parameter \(b=1.5\), using shells spaced \(\Delta r = 1\) apart out to \(r = 3\times b\) gives:

[10]:

from galpy.potential import PlummerPotential, evaluatePotentials
pp= PlummerPotential(amp=1.,b=1.5)
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.],
                             xlabel=r'$r$',ylabel=r'$\Phi(r)$',
                             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=18.,loc='lower right',frameon=False);

../_images/chapters_I-01.-Potential-Theory-and-Spherical-Mass-Distributions_18_0.png

The orange curve shows the actual potential for a Plummer profile (given in Section 3.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\):

[11]:

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.],
                             xlabel=r'$r$',ylabel=r'$\Phi(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=18.,loc='lower right',frameon=False);

../_images/chapters_I-01.-Potential-Theory-and-Spherical-Mass-Distributions_20_0.png

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:

[12]:

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)-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 and using shells out to \(r = 5b\) then gives:

[13]:

rmin, rmax, dr= 0.5,7.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.],
                             xlabel=r'$r$',ylabel=r'$\Phi(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=18.,loc='lower right',frameon=False);

../_images/chapters_I-01.-Potential-Theory-and-Spherical-Mass-Distributions_24_0.png

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 (\(\ref{eq-spherpot}\)) for the potential. Because \(g_r(r) = -\partial \Phi(r) / \partial 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}\,, \end{equation}

or

\begin{equation}\label{eq-spherpot-encmassdens} \Phi(r) = -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 \(\eqref{eq-spherpot-encmassdens}\) is equal to Equation \(\eqref{eq-spherpot}\) 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.

3.3. Circular velocity and dynamical time¶

Next we introduce some of the most basic properties of mass distributions, those of circular velocity and dynamical time. We define these here based on the properties of circular orbits and they therefore require equations for the motion of a body under the influence of gravity. We discuss these equations in far more detail in Chapters 4.1 and 5. All we need for the purposes of this discussion is that for a circular orbit, the centripetal acceleration \(a_r=-v^2/r\) is balanced by the gravitational field. Thus, for a spherical mass distribution using Equation \(\eqref{eq-sphere-gr}\) we have

\begin{equation}\label{eq-centr-balance} a_r = -\frac{v^2}{r} = g_r(r) = -\frac{G\,M(<r)}{r^2}\,, \end{equation}

or

\begin{equation}\label{eq-vcirc} v^2= -r\,g_r(r) = \frac{G\,M(<r)}{r}\,. \end{equation}

This velocity is the circular velocity \(v_c\) at radius \(r\). This equation directly shows that the circular velocity at radius \(r\) is a direct measure of the total mass contained within that radius and measuring the circular velocity as a function of \(r\) therefore directly measures the mass distribution of a spherical distribution. We will see in Chapter 8 that even for very flattened systems, like disk galaxies, the circular velocity is of a similar magnitude as if the flattened mass were distributed spherically. So the circular velocity measures the mass distribution well even for such systems. For reference, the Sun orbits the center of the Milky Way at a distance of \(R_0 \approx 8\,\mathrm{kpc}\) where the circular velocity is \(v_c \approx 220\,\mathrm{km\,s}^{-1}\). Assuming that the mass within the solar orbit is distributed spherically (a rather poor assumption), the values of \(r = 8\,\mathrm{kpc}\) and \(v_c = 220\,\mathrm{km\,s}^{-1}\) mentioned above give an enclosed mass of

\begin{equation}\label{eq-mass-from-vc-at-sun} M(r<8\,\mathrm{kpc}) \approx 9\times 10^{10}\,M_\odot\,. \end{equation}

By comparing to the mass within 8 kpc in a sophisticated mass model for the Milky Way that is shown in Chapter 2.2.3, you can see that this estimate is actually quite accurate!

We can immediately use the relation between the circular velocity and the mass distribution to understand how we know of the presence of dark matter in galaxies (a preview of what we will discuss in detail in Chapter 9). Interstellar gas orbits the centers of galaxies on approximately circular orbits and the velocity of the gas can be measured through Doppler shifts of radio emission lines (like the 21 cm line). This allowed researchers in the 1970s and 1980s to measure the rotation curves—the circular velocity \(v_c(r)\) as a function of \(r\)—at distances \(r\) where few stars and little gas are observed to be present. What they found was that \(v_c\) is almost constant with \(r\). Re-arranging Equation (\(\ref{eq-vcirc}\)) above shows that this implies that the mass distribution behaves as

\begin{equation}\label{eq-mass-vc} M(<r) = \frac{v_c^2\,r}{G} \propto r\,; \end{equation}

the cumulative mass increases linearly with \(r\). Given that little stellar or gaseous material is observed, this implied the existence of a large amount of dark matter in galaxies.

We can define the crossing time or dynamical time as the time necessary to cross the galaxy, \(t_\mathrm{dyn} \approx R /v\). More formally, we may define the dynamical time at \(r\) as the period of a circular orbit at \(r\), the time it takes to orbit the galaxy (or other stellar system) once

\begin{equation}\label{eq-tdyn} t_\mathrm{dyn} = \frac{2\pi\,r}{v_c}\,. \end{equation}

Note that this is not the only possible definition, but all other definitions would involve \(r/v_c\), perhaps with a different pre-factor (see the discussion of the homogeneous sphere below for more justification for this). Using the relation between \(v_c\) and \(M(<r)\) we can express the dynamical time in terms of the average density

\begin{align} t_\mathrm{dyn} & = 2\pi\,r\,\sqrt{\frac{r}{G\,M(<r)}}\nonumber\\ & = \sqrt{\frac{3\pi}{G\,\bar{\rho}}}\,,\label{eq-tdyn-dens} \end{align}

where \(\bar{\rho}\) is the average density within radius \(r\): \(\bar{\rho} = M(<r)/(4\pi r^3/3)\). Thus we have that \(t_\mathrm{dyn} \propto (G\bar{\rho})^{-1/2}\). This formula provides a simple and general way to estimate the dynamical time of any mass distribution, spherical or otherwise. Systems of high average density (such as the centers of galaxies or the centers of globular clusters) have short dynamical times, while low average density regions (such as the outskirts of the halos of galaxies) have long dynamical time scales.

For galaxies, the dynamical times are typically hundreds of Myr to Gyr. For example, using the distance and velocity of the Sun within the Milky Way quoted above, we get that the dynamical time is

\begin{equation} t_{\mathrm{dyn},0} = {2\pi\times 8\,\mathrm{kpc} \over 220\,\mathrm{km\,s}^{-1}} \approx 225\,\mathrm{Myr}\,. \end{equation}

Because galactic rotation curves are close to flat, the dynamical time scales approximately as the distance to the center and, thus, it drops to tens of Myr in the inner kpc and increases to 1 Gyr at \(r \gtrsim 50\) kpc. Galaxies are about 10 Gyr old today, which depending on where you are in the galaxy corresponds to a few dynamical times (in the outskirts) to a few thousand dynamical times (in the innermost regions, with the solar neighborhood being about 40 dynamical times old. Compared to, for example, the \(\approx 4.6\times 10^9\) times that the Earth has gone around the Sun, galaxies are dynamically young.