You are looking at an old version of these notes, please go here for the latest stable version of these notes.

3. Spherical mass distributions

In this chapter we discuss the properties and orbits in spherical mass distributions, starting with some general results from potential theory that apply to any mass distribution. Real galaxies look like this:

Image of M101, a spiral-disk galaxy Image of NGC 660, a polar-ring galaxy

(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. Potential theory

3.1.1. Generalities

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, as we saw in this subsection the gravitational force even between large chunks like stars is dominated by distant bodies, which can be approximated as a smooth density, rather than by the more discrete, nearby bodies. We therefore need to determine the gravitational force \(\vec{F}(\vec{x})\) on a body, like a star, with a mass \(m\), due to a smoothly varying matter density \(\rho(\vec{x})\). We can write the gravitational force \(\vec{F}(\vec{x})\) in terms of 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}

In many applications of galactic dynamics the mass \(m\) of a body interacting with the gravitational field does not matter. This is because intertial mass is equal to gravitational mass in applying Newton’s second law: \(m \vec{a} = \vec{F} = m \vec{g}\) and therefore \(\vec{a} = \vec{g}\). We will use the terms “gravitational force” and “gravitational field” interchangeably and will also often write \(\vec{F}\) for \(\vec{g}\). When we say “force”, we typically mean “field”. Similarly, we will also often use the terms “force” and “acceleration” interchangeably, because \(\vec{a} = \vec{g}\) through Newton’s second law. Obviously we can only do this when the force is that due to gravity.

Sec. 2.1 of Binney & Tremaine (2008; BT08 hereafter) explains that the gravitational field can be written as the gradient of a scalar function \(\Phi(\vec{x})\), the gravitational potential

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

and that the potential is related to the density through the Poisson equation

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

where \(G\) is the gravitational constant. The Poisson equation is the fundamental equation 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”). Because only the derivative of the potential has physical significance (as the gravitational field), we can add or subtract any constant from the potential without changing the dynamics; whenever possible we fix this constant such that the potential equals zero at \(r = \infty\).

3.1.2. Spherical systems: Newton’s theorems

For spherical mass distributions, Newton proved two fundamental theorems that significantly simplify all work with spherical mass distributions. These are:

Newton’s first theorem: A body that is inside a spherical shell of matter experiences no net gravitational force from that shell.

Newton’s second 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.

(See BT08 for proofs of these theorems).

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 6 that this is not the case for flattened mass distributions (the proof of Newton’s first theorem immediately makes it clear why it does not hold for flattened distributions). The first theorem also 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}

The second 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}

This gravitational potential is clearly continuous at \(r=R\). 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 equation for the Laplacian in spherical coordinates from Chapter 2. Thus, for spherical mass distributions, we can compute the potential \(\Phi\) using a simple quadrature for any mass distribution \(\rho(r)\). The only non-zero component of the gravitational field is the radial component \(g_r\), which is given by

\begin{equation} 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}


\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}\) up to a constant using integration by parts. These forms can make calculation of the potential easier if the enclosed mass is known.

3.1.3. Circular velocity, dynamical time, escape velocity

Now consider a body with a negligible mass \(m\)—a test particle—on a circular orbit with radius \(r\) in a spherical mass distribution. The motion of this body is governed by Newton’s second law of motion, where the acceleration is the centripetal acceleration \(a_r=-v^2/r\):

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


\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 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 6 that even for very flattened systems, like disk galaxies, the circular velocity is almost the same 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}\) (e.g., Bovy et al. 2012).

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. Interstellar gas orbits the center 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. 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}

In Chapter 2 we already introduced the crossing time or dynamical time as the time necessary to cross the galaxy, \(t_\mathrm{dyn} \approx R /v\) in the notation of Chapter 2. 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}}\nonumber\\ & = \sqrt{\frac{3\pi}{G\,\bar{\rho}}}\,,\label{eq-tdyn-dens} \end{align}

where \(\bar{\rho}\) is the average density: \(\bar{\rho} = M/(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 a mass distribution. Systems of high density (such as the centers of galaxies or the centers of globular clusters) have short dynamical times, while low density regions (such as the outskirts of the halos of galaxies) have long dynamical time scales.

In most potentials, bodies that move fast enough can “escape” the gravitational field of a mass distribution. This means that such bodies do not orbit the center-of-mass of this mass distribution—because gravity is a long-range force, the gravitational field is technically felt even at extremely large distances. Because the gravitational field can be expressed as the (negative) gradient of a potential, the gravitational force is conservative: the amount of work done against gravity to move a body from position \(\vec{x}_1\) to \(\vec{x}_2\) is independent of the path taken. The energy \(E\) per unit mass of a body in a gravitational field is given by

\begin{equation}\label{eq-energy} E = \frac{|\vec{v}|^2}{2} + \Phi(\vec{x})\,. \end{equation}

For a static potential—a potential that does not change in time—the energy is conserved. This can be shown simply by taking the time derivative of the energy and using Newton’s second law to replace the gradient of the potential with minus the acceleration \(-\dot{\vec{v}}\)

\begin{align} \dot{E} & = \vec{v}\cdot\dot{\vec{v}} + \nabla \Phi(\vec{x})\,\dot{\vec{x}}+\frac{\partial \Phi(\vec{x})}{\partial t}\,\nonumber\\ & = \vec{v}\cdot\dot{\vec{v}}-\dot{\vec{v}}\cdot\vec{v} +\frac{\partial \Phi(\vec{x})}{\partial t}\,\label{eq-energy-conserved}\\ & = \frac{\partial \Phi(\vec{x})}{\partial t}\,,\nonumber\\ \end{align}

which is zero for a static potential (a static potential does not explicitly depend on time and therefore \(\partial \Phi / \partial t = 0\)). That the energy is conserved allows us to compute the escape velocity \(v_\mathrm{esc}\): the minimum velocity necessary to “escape” the gravitational field, starting from position \(r\) to a resting position at \(r=\infty\) (assuming a spherical mass distribution). Energy conservation gives

\begin{equation}\label{eq-vesc} \frac{v_\mathrm{esc}^2}{2} + \Phi(r) = \frac{[v=0]^2}{2} + \Phi(\infty)\,, \end{equation}


\begin{equation} v_\mathrm{esc} = \sqrt{2[\Phi(\infty)-\Phi(r)]}\,, \end{equation}

Because we normally try to define the potential such that \(\Phi(\infty)=0\), this typically simplifies to \(v_\mathrm{esc} = \sqrt{-2\Phi(r)}\). Thus, the escape velocity measures the depth of the potential well. This means that the escape velocity contains information about the mass distribution outside of the radius at which it is measured, unlike the circular velocity, which only depends on the mass contained within the radius at which it is measured. This should be clear from the meaning of the “escape” velocity: the speed necessary for a body to escape must depend on the cumulative amount of gravitational attraction the body needs to overcome out to \(\infty\). It is also clear from the expression for the gravitational potential in Equation (\(\ref{eq-spherpot}\)): the second term depends on the density outside of \(r\).

In the Milky Way it is possible to measure the escape velocity near the Sun by looking at the distribution of velocities and attempting to identify a high-velocity cut-off. The cut-off arises because stars at the escape velocity and beyond should be very rare, as they will escape the Galaxy and never return. This has been done with ever-increasing data from large-scale surveys measuring velocities for solar neighborhood stars (e.g., Leonard & Tremaine 1990; Smith et al. 2007; Piffl et al. 2014) with a local escape speed of \(v_\mathrm{esc} \approx 550\,\mathrm{km\,s}^{-1}\). We can get a rough estimate of the total mass of the Milky Way by comparing this measured escape velocity to the escape velocity that the Milky Way would have if all of its mass were within the solar orbit. We can use Newton’s second theorem and Equation (\(\ref{eq-spherpot-outside}\)) to get the escape velocity

\begin{equation} v_\mathrm{esc}(r=8\,\mathrm{kpc}) = \sqrt{2\frac{GM(r<8\,\mathrm{kpc})}{8\,\mathrm{kpc}}}\,, \end{equation}

and substitute the estimate of \(M(r<8\,\mathrm{kpc})\) from Equation (\(\ref{eq-mass-from-vc-at-sun}\)). This gives \(v_\mathrm{esc} = 311\,\mathrm{km\,s}^{-1}\), well below the observed value of \(550\,\mathrm{km\,s}^{-1}\). Equation (\(\ref{eq-spherpot}\)) demonstrates that we need to assume a density law outside of the solar orbit if we want to match the observed escape velocity (the potential outside of \(r\) depends on \(\int_r^\infty\mathrm{d}r'\,\rho(r')\,r'\), not the mass). Let’s use the mass–radius relation that we derived from the fact that rotation curves are flat in Equation (\(\ref{eq-mass-vc}\)): \(\rho \propto r^{-2}\). For this particular density law, the mass and potential do not converge as we go to \(r\rightarrow \infty\), so let’s define the edge of the Milky Way as being at \(100\,\mathrm{kpc}\) outside of which the density is zero (this is just a rough guess, given that we need to add lots of mass to explain the measured escape velocity; later we will see that the standard definition of “edge” gives a radius around \(250\,\mathrm{kpc}\)). Then we can show that the potential difference between \(r_\infty = 100\,\mathrm{kpc}\) and \(r_0 = 8\,\mathrm{kpc}\) is given in terms of the mass \(M(r<8\,\mathrm{kpc})\) and the mass \(\Delta M\) between \(r_0\) and \(r_\infty\)

\begin{align} \Delta \Phi & = \Phi(r_\infty)-\Phi(r_0)\\ & = -G\,\left[\frac{M(r<8\,\mathrm{kpc})+\Delta M}{r_\infty}-\frac{M(r<8\,\mathrm{kpc})}{r_0}-\frac{\Delta M}{r_\infty-r_0}\,\ln (r_\infty/r_0)\right]\nonumber\,. \end{align}

Plugging this into \(v_\mathrm{esc} = \sqrt{2\Delta\Phi}\) and using the enclosed mass from Equation (\(\ref{eq-mass-from-vc-at-sun}\)) we find that (this calculation can be simplified by assuming that \(M(r<8\,\mathrm{kpc}) \ll \Delta M\))

\begin{equation} \Delta M \approx 1.4\times 10^{12}\,M_\odot\,! \end{equation}

That the escape velocity near the Sun is approximately \(550\,\mathrm{km\,s}^{-1}\) therefore implies that the total mass of the Milky Way must be \(\approx 10^{12}\,M_\odot\), more than an order of magnitude larger than the mass we estimate to lie within the solar orbit. This simple estimate is remarkably accurate: a variety of complicated methods show that the Milky Way’s total mass is \(\approx 0.8\) to \(1.6 \times 10^{12}\,M_\odot\), in good agreement with our estimate. This mass is again much larger than the total mass that we see in stars and gas, again directly implying the presence of a large amount of dark matter in the Milky Way.

3.2. Examples of spherical potentials

Let’s now consider some example spherical potentials that illustrate the concepts above. This is by no means an exhaustive list of spherical potentials. Sec. 2.2.2 of BT08 list a few more examples or look at the spherical potentials included in galpy for some more examples.

3.2.1. Point-mass potential

The most basic spherical potential is that of a point mass. Because every point is outside of the (zero radius) “shell” defined by the point mass, Equation \(\eqref{eq-spherpot-outside}\) tells us that the potential is

\begin{equation}\label{eq-pot-kepler} \Phi(r) = -\frac{G\,M}{r}\,, \end{equation}

for a point mass \(M\). The circular velocity is given by

\begin{equation}\label{eq-vc-kepler} v_c(r) = \sqrt{\frac{G\,M}{r}}\,. \end{equation}

The escape velocity is left for an “end-of-chapter” question. The point-mass potential of the Sun is the dominant gravitational potential that sets planetary trajectories in the solar system (perturbations from other planets are very small). Because Equation \(\eqref{eq-vc-kepler}\) is equivalent to Kepler’s third law, first proposed by Kepler, the point-mass potential and orbits within it are often referred to as Keplerian. For example, a disk of mass \(m \ll M\) orbiting a point mass is “Keplerian” and is in “Keplerian rotation”. In galactic dynamics, point-mass potentials are typically only used to represent the gravitational field from supermassive black holes or as a crude approximation. Point-mass potentials are not typically used in \(N\)-body simulations of galaxies as we will discuss in Chapter 10: the gravitational potential of individual particles in \(N\)-body simulations is smoothed to, for example, a Plummer potential (see below) to avoid (unphysical) close encounters; these encounters are unphysical because galaxies are collisionless (see Chapter 2.3). \(N\)-body simulations of stellar clusters do use the point-mass potentials for all of their constituent stars.

3.2.2. The homogeneous sphere

The potential of a homogeneous density sphere, \(\rho(r) = \rho_0 =\) constant can be calculated using \(\eqref{eq-spherpot}\); the result is

\begin{equation}\label{eq-pot-homogeneous} \Phi(r) = \frac{2\pi\,G\,\rho_0}{3}\,r^2\,, \end{equation}

where we have dropped a constant term equal to \(\infty\). This quadratic potential is that of a harmonic oscillator. The circular velocity is

\begin{equation} v_c(r) = \sqrt{\frac{4\pi\,G\,\rho_0}{3}}\,r\,, \end{equation}

and the dynamical time using the definition from \(\eqref{eq-tdyn}\) is therefore

\begin{equation}\label{eq-tdyn-homogeneous} t_\mathrm{dyn} = \sqrt{\frac{3\pi}{G\,\rho_0}}\,. \end{equation}

The dynamical time in Equation \(\eqref{eq-tdyn-homogeneous}\) is the same as that we derived in Equation \(\eqref{eq-tdyn-dens}\) when substituting \(\rho_0\) for the average density. In a homogeneous sphere, the dynamical time is constant with \(r\). If instead of considering a body on a circular orbit, we look at a body moving perfectly radially, it’s equation of motion according to Newton’s second law is

\begin{equation} \ddot{r} = -\frac{4\pi\,G\,\rho_0}{3}\,r\,, \end{equation}

which is the equation of motion of a simple harmonic oscillator with frequency \(\omega = \sqrt{4\pi\,G\,\rho_0/3}\). The period of this radial oscillation is

\begin{equation} T = 2\pi/\omega \propto (G\,\rho_0)^{-1/2} \propto t_\mathrm{dyn} = \mathrm{constant}\,, \end{equation}

which is again constant. This result provides further justification for the definition of the dynamical time in \(\eqref{eq-tdyn}\): no matter whether a body is on a circular orbit or a purely radial orbit, the period of the orbit is \(\approx t_\mathrm{dyn}\).

Because \(v_c(r) \propto r\) for a homogeneous sphere, a disk of bodies rotates with an angular rotation rate of \(\Omega(r) = v_c(r)/r =\) constant. Because the rotation rate is constant, such a disk rotates as a solid body and this setup is referred to as solid-body rotation.

3.2.3. Plummer sphere

Next we will look at a few generalizations of the Kepler point-mass potential above. One of the most important of these is the Plummer model, which smooths the Kepler potential by replacing the radius \(r\) in the denominator of Equation \(\eqref{eq-pot-kepler}\) by \(\sqrt{r^2+b^2}\) with \(b\) a constant, the Plummer scale length,

\begin{equation}\label{eq-pot-plummer} \Phi(r) = -\frac{G\,M}{\sqrt{r^2+b^2}}\,. \end{equation}

This potential approaches a Kepler potential as \(b \rightarrow 0\), or, more usefully, \(r \gg b\). This potential was originally introduced to represent globular clusters, but is no longer a preferred models for these systems. It remains highly useful because it has many convenient analytical properties and provides a crude model for, e.g., dark matter halos that is easy to work with and it is also used to represent the stellar density in small dwarf galaxies. It is also a simple kernel used to smooth the gravitational field of point particles in \(N\)-body simulations.

It is straightforward to derive the circular velocity, escape velocity, and density profile from Equation \(\eqref{eq-pot-plummer}\) (using the Poisson equation to obtain the density). See BT08 for further details.

3.2.4. Isochrone potential

Similar to the Plummer sphere, the isochrone potential is obtained by smoothing the potential of a point mass. For the isochrone potential, this is done by replacing \(r\rightarrow b+\sqrt{r^2+b^2}\) in the denominator of the gravitational potential (see Equation \([\ref{eq-pot-kepler}]\)), with \(b\) again a constant

\begin{equation}\label{eq-pot-isochrone} \Phi(r) = -\frac{G\,M}{b+\sqrt{r^2+b^2}}\,. \end{equation}

Similar to the Plummer sphere, this potential approaches a Kepler potential as \(b \rightarrow 0\), or as \(r \gg b\). We can determine how fast the isochrone model approaches the Kepler potential by expaning the potential around \(b/r\) for \(r \gg b\)

\begin{align} \Phi_{r\gg b}(r) & = -\frac{G\,M}{r\left[b/r+\sqrt{1+\left(b/r\right)^2}\right]}\nonumber\\ & \approx -\frac{G\,M}{r\left[b/r+1+\frac{1}{2}\,\left(b/r\right)^2\right]}\\ & \approx -\frac{G\,M}{r\left[1+b/r\right]}\nonumber\\ & \approx -\frac{G\,M}{r}\,\left(1-\frac{b}{r}\right)\,.\nonumber\\ \end{align}

Because the first term is the potential of a point mass, which has zero density at \(r\gg b\), the density at large radii is given by that corresponding to the second term, \(\Phi(r) \propto r^{-2}\). From the Poisson equation, we know that the density is a second derivative of the potential, and must therefore tend to \(\rho_{r\gg b}(r) \propto r^{-4}\). This \(1/r^4\) behavior is also that of a popular model for dark-matter halos, the Hernquist model. We will discuss power-law potentials and the Hernquist model in more detail below.

When \(r/b \rightarrow 0\), we can expand the potential around \(r/b\)

\begin{align} \Phi_{r\ll b}(r) & = -\frac{G\,M}{b\left[1+\sqrt{\left(r/b\right)^2+1}\right]}\nonumber\\ & \approx -\frac{G\,M}{b\left[1+1+\frac{1}{2}\,\left(r/b\right)^2\right]}\\ & \approx -\frac{G\,M}{2b}\,\left(1-\frac{1}{4}\,\frac{r^2}{b^2}\right)\nonumber\\ & = \frac{G\,M}{8\,b^3}\,r^2+\mathrm{constant}\nonumber\,. \end{align}

Up to an irrelevant constant, this is the potential of a homogeneous sphere with density \(\rho_0 = 3\,M/(16\pi\,b^3)\). This immediately tells us that the central density of the isochrone mass distribution is this \(\rho_0\).

Thus, the isochrone potential smoothly interpolates between the potential of a point mass and that of a homogenous density distribution, essentially the two extremes of the types of densities seen in stellar systems.

The circular velocity of the isochrone potential is given by

\begin{equation}\label{eq-vc-isochrone} v_c^2(r) = \frac{G\,M\,r^2}{(b+a)^2\,a}\,, \end{equation}

where \(a = \sqrt{r^2+b^2}\).

It is clear that the isochrone model can represent a wide range of potentials. Because it interpolates between the two extremes of density distributions, over a narrow radial range it can essentially represent all spherical potentials. But the main reason of the importance of the isochrone potential is that it is the most realistic model for galaxies in which all of the orbits can be solved analytically. That is, just like for the Kepler potential, we can work out the full orbit of a body in terms of elementary functions, without requiring any numerical quadrature or numerically solving the differential equation of motion. Naturally, the analytic solution is more complicated than that for a Kepler potential, but not so much more so. That all of the orbits are analytic is extremely useful in many settings in galactic dynamics.

Let’s compare the rotation curves of the spherical potentials that we have discussed so far by plotting them with galpy. We normalize the point-mass, Plummer, and isochrone potentials such that they all have \(G\,M = 1\), and the homogeneous sphere such that \(v_c = 1\) for \(r = 1\) (the latter, because the homogeneous sphere does not have finite mass). Because galpy does not actually contain a pure homogeneous sphere potential, we approximate it by using an isochrone potential with \(b \gg 1\). We set the scale \(b\) of the isochrone and Plummer potentials to \(b=1\).

In [5]:
from galpy import potential
from galpy.util import bovy_plot
# Define each potential
kp= potential.KeplerPotential(amp=1.)
hp= potential.IsochronePotential(normalize=1.,b=1000.)
pp= potential.PlummerPotential(amp=1.,b=1.)
ip= potential.IsochronePotential(amp=1.,b=1.)
# Plot the rotation curve for each
line_kepler= potential.plotRotcurve(kp,Rrange=[0.,10.],\
                label=r'$\mathrm{Point\ mass}$',yrange=[0.,1.5],
line_homog= potential.plotRotcurve(hp,Rrange=[0.,10.],\
                label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
line_plummer= potential.plotRotcurve(pp,Rrange=[0.,10.],\
line_isochrone= potential.plotRotcurve(ip,Rrange=[0.,10.],\
       fontsize=18.,loc='upper right',frameon=False);

We see the expected behavior in these curves: the point-mass \(v_c(r)\) drops as \(r^{-1/2}\), while the other potentials all have \(v_c(r) \propto r\) at \(r \ll b\) (and everywhere for the homogeneous sphere, which quickly runs out of the shown range in \(v_c\)). At \(r > b\) the Plummer and isochrone potentials approach the Keplerian curve, because they all have the same total mass. The Plummer potential approaches the Keplerian curve quicker than the isochrone potential for the same \(b\).

3.2.5. Power-law models

An important class of models are those in which the density is a power-law of radius

\begin{equation} \rho(r) = \rho_0\,\left(\frac{r_0}{r}\right)^\alpha\,, \end{equation}

where \(\rho_0\) is the density at \(r=r_0\). The enclosed mass of this density profile is

\begin{equation} M(<r) = \frac{4\pi\,\rho_0\,r_0^3}{3-\alpha}\,\left(\frac{r_0}{r}\right)^{\alpha-3}\,, \end{equation}

which is finite at \(r=\infty\) only when \(\alpha < 3\). Using Equation \(\eqref{eq-spherpot-encmass}\), we find that the potential is

\begin{equation} \Phi(r) = -\frac{4\pi\,G\,\rho_0\,r_0^2}{(3-\alpha)(2-\alpha)}\,\left(\frac{r_0}{r}\right)^{\alpha-2}\,,\qquad (\alpha \neq 2)\,, \end{equation}


\begin{equation}\label{eq-pot-powlaw-log} \Phi(r) = \frac{4\pi\,G\,\rho_0\,r_0^2}{(3-\alpha)}\ln r\,,\qquad (\alpha = 2)\,. \end{equation}

The circular velocity is given by

\begin{equation} v^2_c(r) = \frac{4\pi\,G\,\rho_0\,r_0^2}{3-\alpha}\,\left(\frac{r_0}{r}\right)^{\alpha-2}\,. \end{equation}

For \(\alpha = 2\) the circular velocity is constant \(v_c\) and we can write the gravitational potential in Equation \(\eqref{eq-pot-powlaw-log}\) in terms of \(v_c\)

\begin{equation}\label{eq-pot-powlaw-log-vc} \Phi(r) = v_c^2\,\ln r\,,\qquad (\alpha = 2)\,. \end{equation}

This is an important potential that is, for obvious reasons, often referred to as the logarithmic potential. Because it gives rise to a constant circular velocity—a flat rotation curve—it is a very simple and often-used model for the potential of galaxies with a flat rotation curve. When we study self-consistent dynamical equilibrium states in Chapter 4, we will see that for the logarithmic potential the velocity dispersion is constant with \(r\) and it is therefore also often referred to as the singular isothermal sphere, with the “singular” modifier due to the \(1/r^2\) density divergence as \(r\rightarrow 0\).

3.2.6. Two-power density models

A final important class of spherical potentials for galaxies are those that are described by a different power law in their inner and outer parts, with a smooth transition region. A popular set of models with this property is that with the following density law

\begin{equation}\label{eq-twopower-dens} \rho(r) = \frac{\rho_0\,a^\alpha}{r^\alpha\,(1+r/a)^{\beta-\alpha}}\,. \end{equation}

In the small and large \(r\) limit, these models behave as power law models

\begin{equation} \rho(r) = \begin{cases} \rho_0\,\left(\frac{a}{r}\right)^\alpha, & r \ll a\,,\\ \rho_0\,\left(\frac{a}{r}\right)^\beta, & r \gg a\,. \end{cases} \end{equation}

The inner slope is therefore \(\alpha\) and the outer slope is \(\beta\). We will only consider two popular instances of this set of models: the Hernquist model, which has \(\alpha = 1, \beta =4\) and the Navarro-Frenk-White (NFW) model, with \(\alpha = 1, \beta = 3\). See Sec. 2.2.2(g) of BT08, Dehnen (1993), and Tremaine et al. (1994) for more information on this set of models. The enclosed mass is given by

\begin{equation}\label{eq-massenc-hernnfw} M(<r) = 4\pi\,\rho_0\,a^3\,\times \begin{cases} \frac{(r/a)^2}{2\,(1+r/a)^2}, & \mathrm{(Hernquist)}\,\\ \ln(1+r/a)-\frac{r/a}{1+r/a}, & \mathrm{(NFW)}\,. \end{cases} \end{equation}

Using Equation \(\eqref{eq-spherpot-encmass}\), this gives the gravitational potential

\begin{equation}\label{eq-potential-hernnfw} \Phi(r) = -4\pi\,G\,\rho_0\,a^2\,\times \begin{cases} \frac{1}{2\,(1+r/a)}, & \mathrm{(Hernquist)}\,\\ a\,\ln(1+r/a)/r, & \mathrm{(NFW)}\,. \end{cases} \end{equation}

The total mass of an NFW potential does not converge as we increase \(r\), but Equation \(\eqref{eq-massenc-hernnfw}\) demonstrates that the total mass of a Hernquist potential is \(M = 2\pi\,\rho_0\,a^3\). The Hernquist potential written in terms of this mass \(M\) is

\begin{equation} \Phi(r) = -\frac{G\,M}{r+a}\,. \end{equation}

The Hernquist potential is therefore also a simple smoothing of the Kepler point-mass potential (like the Plummer model), in this case by replacing \(r \rightarrow r+a\). As for the Plummer model, this often means that calculations involving a Hernquist potential can be analytically simplified.

The NFW profile (Navarro et al. 1997) is found to describe the density profile of dark matter halos well in cosmological simulations (Navarro et al. 1996) and remains highly important. Hernquist profiles have somewhat more tractable properties and are also a popular model for dark matter halos and for elliptical galaxies and galactic bulges, because it can approximately represent the \(r^{1/4}\) de Vaucouleurs profile (see Chapter 2). The outskirts of dark-matter halos at the present time have not yet fully formed and it has also been argued that the long-term state of dark-matter halos more closely resembles a Hernquist profile than an NFW profile (Busha et al. 2005).

Let’s compare the rotation curves of Hernquist and NFW potentials with the same scale radius \(a\). First, we compare the rotation curves for the case where both potentials correspond to the same enclosed mass at \(r= 12\,a\) (this is approximately the virial radius for Milky-Way-mass dark-matter halos). We do this by computing the enclosed mass using galpy’s .mass function for each potential set up with an amplitude of unity and then setting the new amplitude such that both potentials have the same mass:

In [6]:
from galpy import potential
from galpy.util import bovy_plot
# Define each potential, equal mass at 12 scale radii
hp_4scale= potential.HernquistPotential(amp=1.,a=1.)
hp= potential.HernquistPotential(amp=10./hp_4scale.mass(12.),a=1.)
nfp_4scale= potential.NFWPotential(amp=1.,a=1.)
nfp= potential.NFWPotential(amp=10./nfp_4scale.mass(12.),a=1.)

# Plot the rotation curve for each
line_hernquist= potential.plotRotcurve(hp,Rrange=[0.,15.],\
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
       fontsize=18.,loc='upper right',frameon=False)
bovy_plot.bovy_text(r'$\mathrm{Equal\ mass\ at}\ R=12\,a$',

Because both potentials are defined to have the same mass at \(r=12\,a\), the circular velocity at these radii is the same. We see that the rotation curve of the Hernquist potential reaches a much higher maximum value than the NFW potential with the same mass. At \(r > 12\,a\), the Hernquist \(v_c(r)\) declines more steeply with radius than the NFW \(v_c(r)\), because of the \(r^{-4}\) behavior of the density at \(r \gg a\) for the Hernquist profile versus the \(r^{-3}\) behavior for the NFW density.

Next, let’s compare the rotation curves if we normalize the Hernquist and NFW potentials such that they have the same inner density profile \(\rho(r) \propto \rho_0\,r^{-1}\) at \(r \ll a\):

In [7]:
from galpy import potential
from galpy.util import bovy_plot
# Define each potential, equal inner profile
hp= potential.HernquistPotential(amp=20.,a=1.)
nfp= potential.NFWPotential(amp=20.,a=1.)
# Plot the rotation curve for each
line_hernquist= potential.plotRotcurve(hp,Rrange=[0.,15.],\
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
       fontsize=18.,loc='upper right',frameon=False)
bovy_plot.bovy_text(r'$\mathrm{Equal\ inner\ density\ profile}$',

The inner rotation curves are now the same, because the enclosed mass profile is the same at \(r \ll a\), but because the density of the NFW profile is shallower than that of the Hernquist profile, the NFW rotation curve reaches a significantly higher value than that of the Hernquist profile.

3.3. Orbits in spherical potentials

“Orbits” are the trajectories that bodies travel on under the influence of gravity in a gravitational field. These trajectories can be computed using Newton’s second law

\begin{equation}\label{eq-newton-eom} m\, \ddot{\vec{x}} = m\,\vec{g}(\vec{x})\,, \end{equation}

where \(\vec{g}(\vec{x})\) is the gravitational field, \(m\) is the mass of the body, and the double dot \(\ddot{\vec{x}}\) indicates the second time derivative of the three-dimensional position \(\vec{x}\). As discussed above, because we can cancel the mass on both sides of this equation, the acceleration \(\ddot{\vec{x}}\) is independent of the mass of the body (which is not the case, for example, for the acceleration due to the electromagnetic force—the Lorentz force). Because of this property, we typically ignore the actual mass \(m\) of bodies orbiting on a gravitational field and consider familiar quantities such as the energy, momentum, angular momentum, and the force per unit mass. For example, the momentum per unit mass is \(\dot{\vec{x}}\) and we will typically refer to this simply as the “momentum”.

3.3.1. General properties of orbits in spherical potentials

Equation \(\eqref{eq-newton-eom}\) is a ordinary differential equation that typically needs to be solved numerically. Only for a handful of galactic potentials can Equation \(\eqref{eq-newton-eom}\) be solved analytically. Even for most spherical potentials, the solution of Equation \(\eqref{eq-newton-eom}\) needs to be numerical.

Sec. 3.1 of BT08 describes general results for orbits in spherical potentials. We summarize the main results here, but read Sec. 3.1 in its entirey (except for 3.1(d) hyperbolic encounters) to get the full picture. Because we are dealing with spherical mass distributions, we will mostly work in spherical coordinates and will denote the three-dimensional position as \(\vec{r}\) rather than \(\vec{x}\). Because only the radial derivative of the potential is non-zero, the equation of motion in a spherical potential simplifies to

\begin{equation}\label{eq-eom-spher} \ddot{\vec{r}} = g_r(r)\,\hat{\vec{r}}\,, \end{equation}

where \(\hat{\vec{r}}\) is the radial unit vector and \(\vec{r} = r\,\hat{\vec{r}}\) (with \(r = |\vec{r}|\) the magnitude of the radius). The force points along the radial direction.

Spherical symmetry implies invariance under any rotation around the center of the potential. By Noether’s theorem this implies the existence of a conserved quantity. This quantity is the angular momentum vector

\begin{equation} \vec{L} = \vec{r}\times \dot{\vec{r}}\,, \end{equation}

the cross product of the position vector and the velocity vector. That the angular momentum vector is conserved is most easily demonstrated simply by taking its derivative

\begin{align} \dot{\vec{L}} & = \frac{\mathrm{d}}{\mathrm{d} t}\left(\vec{r}\times\dot{\vec{r}}\right)\nonumber\\ & = \dot{\vec{r}}\times\dot{\vec{r}} + \vec{r}\times\ddot{\vec{r}}\\ & = g_r(r)\,r\,\,\hat{\vec{r}}\times\hat{\vec{r}} = 0\nonumber\,, \end{align}

where the first term in the second line is zero because the cross product of any vector with itself is zero. The third line uses Equation \(\eqref{eq-eom-spher}\) and \(\vec{r} = r\,\hat{\vec{r}}\); the cross product in the third line is therefore again that of a vector with itself. That the angular momentum vector is conserved implies that the orbit is constrained to a single plane and we therefore only have to consider the motion in this plane. The magnitude of the angular momentum is \(L\). Using polar coordinates \((r,\psi)\) in this plane, the relevant equations of motion in this plane are

\begin{align}\label{eq-eom-spher-r} \ddot{r} = \frac{L^2}{r^3} + & g_r(r) = \frac{L^2}{r^3} - \frac{\mathrm{d} \Phi(r)}{\mathrm{d} r}\,\\ \dot{\psi} & = \frac{L}{r^2}\,.\label{eq-eom-spher-angle} \end{align}

The energy \(E\) (see Equation [\(\ref{eq-energy}\)]) of an orbit in a spherical potential is given by

\begin{align} E & = \frac{\dot{r}^2}{2} + \frac{r^2\,\dot{\psi}^2}{2} + \Phi(r)\nonumber\\ & = \frac{\dot{r}^2}{2} + \frac{L^2}{2\,r^2} + \Phi(r)\,\label{eq-energy-spherpot} \end{align}

where we have replaced \(\dot{\psi}\) with \(L\) using Equation \(\eqref{eq-eom-spher-angle}\). As shown in Equation \(\eqref{eq-energy-conserved}\), this energy is conserved if the potential \(\Phi(r)\) does not explicitly depend on time. The equation of motion for \(r\) in Equation \(\eqref{eq-eom-spher-r}\) can be re-written in terms of an effective potential \(\Phi_\mathrm{eff}(r)\) given by

\begin{equation} \Phi_\mathrm{eff}(r) = \Phi(r) + \frac{L^2}{2\,r^2}\,, \end{equation}

and the energy from Equation \(\eqref{eq-energy-spherpot}\) is then simply

\begin{equation} E = \frac{\dot{r}^2}{2} + \Phi_\mathrm{eff}(r)\,. \end{equation}

The energy can therefore be thought of as a radial energy.

The energy and angular momentum work together to restrict the part of physical space that an orbit can occupy. Because for typical smooth galactic potentials \(\Phi(r)\) increases with increasing \(r\), energy conservation restricts orbits with energy \(E\) to those \(r\) with \(\Phi(r) \leq E\), which restricts bound orbits to \(r \leq r_{\mathrm{max}}(E)\) for which \(\Phi(r_\mathrm{max}) \leq E\), but allows orbits to reach \(r=0\); unbound orbits have \(E > \Phi(\infty)\) and therefore have no \(r_{\mathrm{max}}(E)\). The conservation of angular momentum in a spherical potential changes this restriction to \(\Phi_\mathrm{eff}(r) \leq E\). Like the potential itself, the effective potential decreases with decreasing \(r\) at large \(r\), but increases with decreasing \(r\) at small \(r\), because of the \(L^2/[2r^2]\) term, called the angular momentum barrier. Any spherical potential corresponding to a spherical density at most decreases as fast as the Kepler potential, \(\Phi(r) \propto 1/r\). Therefore, the increase in \(\Phi_\mathrm{eff}(r)\) at small \(r\) from the angular-momentum barrier \(L^2/[2r^2]\) always becomes larger than the decrease from \(\Phi(r)\) at small enough \(r\) and the \(L^2/[2r^2]\) term grows without bounds as \(r\rightarrow 0\). Therefore, for bound orbits, the equation \(\Phi_\mathrm{eff}(r) = E\) has two solutions, one at small \(r\) and one at large \(r\). Unbound orbits have no maximum \(r\), but as long as they have non-zero angular momentum, they do have a minimum \(r\).

Bound orbits therefore perform oscillations in \(r\) between a minimum and maximum value. Analogous to the perihelion and apohelion distances of orbits in the Keplerian potential for the solar system, these are known as the pericenter \(r_p\) and the apocenter \(r_a\), respectively. Because these radii are minima and maxima of \(r(t)\), we can find them by solving \(\dot{r} = 0\). Using Equation \(\eqref{eq-energy-spherpot}\), we can easily solve for these by solving \(\dot{r}^2 = 0\)

\begin{equation}\label{eq-spherpot-drdt} \dot{r}^2 = 2[E-\Phi(r)]-\frac{L^2}{r^2} = 0\,, \end{equation}

which is the same as solving \(\Phi_\mathrm{eff}(r) = E\), but makes it clear that these radii are turning points of the orbit. In general, this equation needs to be solved numerically. Circular orbits have a constant radius and therefore \(r_p = r_a\). A convenient quantity to estimate how circular an orbit is is given by the orbital eccentricity \(e\) which is defined as

\begin{equation}\label{eq-ecc-spher} e = \frac{r_a-r_p}{r_a+r_p}\,. \end{equation}

For a circular orbit \(e=0\), while orbits that are close to unbound have \(r_a \gg r_p\) and therefore \(e \rightarrow 1\) as \(r_a\) grows.

The time necessary to travel from pericenter to apocenter and back is known as the radial period \(T_r\). \(T_r\) can be computed as

\begin{align} T_r & = 2\,\int_{r_p}^{r_a}\mathrm{d} t\nonumber\\ & = 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{\mathrm{d} t}{\mathrm{d} r}\\ & = 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{1}{\sqrt{2[E-\Phi(r)]-L^2/r^2}}\nonumber\,,\\ \end{align}

where we used Equation \(\eqref{eq-spherpot-drdt}\). The radial period describes the motion in the radial direction only. The motion of a body in a spherical potential in the azimuthal \(\psi\) direction also has a period associated with it, the azimuthal period \(T_\psi\). This is most easily computed by first calculating the change \(\Delta \psi\) in \(\psi\) over a radial period and then computing the number of radial periods necessary to complete one full \(\psi\) rotation. We have that

\begin{align} \Delta \psi & = 2\,\int_{r_p}^{r_a}\mathrm{d} \psi\nonumber\\ & = 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{\mathrm{d}t}{\mathrm{d}r}\,\frac{\mathrm{d}\psi}{\mathrm{d}t}\\ & = 2\,L\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{1}{r^2\,\sqrt{2[E-\Phi(r)]-L^2/r^2}}\nonumber\,,\\ \end{align}

where we used Equation \(\eqref{eq-spherpot-drdt}\) and Equation \(\eqref{eq-eom-spher-angle}\) to substitute \(\mathrm{d}r/\mathrm{d}t\) and \(\mathrm{d}\psi/\mathrm{d}t\), respectively. The azimuthal period is then

\begin{equation}\label{eq-tr-tpsi-dpsi} T_\psi = \frac{2\pi}{|\Delta \psi|}\,T_r\,, \end{equation}

where we need the absolute value of \(\Delta \psi\), because \(\Delta \psi\) can be negative. All of the quantities \(T_r\), \(\Delta \psi\), and \(T_\psi\) typically need to be computed numerically.

3.3.2. Orbits in specific spherical potentials

To illustrate the concepts introduced in the previous subsection, we consider orbital properties in a few simple spherical potentials. Orbits in the homogeneous sphere

The homogeneous sphere (see Section 3.2.2 above) has particularly simple orbits. For the homogeneous sphere it is simpler to solve the equations of motion in cartesian coordinates rather than in polar coordinates as in Equation \(\eqref{eq-eom-spher-r}\). First, we re-write the potential for the homogeneous sphere (Equation \([\ref{eq-pot-homogeneous}]\)) as

\begin{equation} \Phi(r) = \frac{1}{2}\,\omega^2\,r^2\,, \end{equation}

with \(\omega^2 = 4\pi\,G\,\rho_0/3\). In Cartesian coordinates this becomes

\begin{equation} \Phi(x,y) = \frac{1}{2}\,\omega^2\,\left(x^2+y^2\right)\,. \end{equation}

The equations of motion are therefore

\begin{align} \ddot{x} &= -\omega^2\,x\,,\\ \ddot{y} &= -\omega^2\,y\,. \end{align}

These have the solution

\begin{align} x(t) &= X\,\cos\left(\omega t + \phi_x\right)\,,\\ y(t) &= Y\,\cos\left(\omega t + \phi_y\right)\,. \end{align}

This solution describes an ellipse, with the center at the origin of the mass distribution (\(r=0\)). Note that these are not Keplerian orbits, even though Keplerian orbits are also ellipses; as we will see below, Keplerian orbits are ellipses with the origin of the mass distribution at the focus of the ellipse, not at the center of the ellipse as for the homogeneous sphere. This solution has four constants that need to be determined from the initial condition of the orbit (remember that we are working in the orbital plane, which is constant, and there are therefore only two positions and velocities that count; the position and velocity perpendicular to this plane are zero). These constants are (i+ii) the semi-major and semi-minor axes of the ellipse, \(a\) and \(b\) respectively, (iii) the orientation of the ellipse within the plane, and (iv) the position along the orbit at \(t=0\). Let’s fix the latter two of these by orienting the ellipse such that the semi-major axis is along \(x\) and the orbit is at \(y=0\) at positive \(x\) at \(t=0\):

\begin{align}\label{eq-homogeneous-orbit-1} x(t) & = a\,\cos \omega t\,,\\ y(t) & = b\,\sin \omega t\,.\label{eq-homogeneous-orbit-2} \end{align}

It is then straightforward to demonstrate that the energy and angular momentum are

\begin{align} E & = \frac{\left(a^2+b^2\right)\,\omega^2}{2}\,,\\ L & = a\,b\,\omega\,, \end{align}

or in terms of the eccentricity \(\varepsilon\) of the ellipse: \(\varepsilon = \sqrt{1-b^2/a^2}\)

\begin{align}\label{eq-homogeneous-E} E & = a^2\,\left(1-\frac{\varepsilon^2}{2}\right)\,\omega^2\,,\\ L & = a^2\,\sqrt{1-\varepsilon^2}\,\omega\,.\label{eq-homogeneous-L} \end{align}

The peri- and apocenter radii are simply

\begin{align} r_p & = a\,,\\ r_a & = b\,, \end{align}

which directly follows from the properties of ellipses. Using Equations \(\eqref{eq-homogeneous-E}\) and \(\eqref{eq-homogeneous-L}\) it is also easy to show that \(E = \omega^2 a^2 / 2 + L^2 / [2a^2] = \omega^2 b^2/2 + L^2 / [2b^2]\). Note that the general definition of the orbital eccentricity \(e\) in Equation \(\eqref{eq-ecc-spher}\) does not agree with the standard definition of the eccentricity of an ellipse \(\varepsilon\) for orbits in the homogeneous sphere: \(e \neq \varepsilon\). We will see below that for orbits in the point-mass potential the two definitions to agree.

Because an ellipse goes through its pericenter and apocenter twice for each full rotation around the center, the radial period is half of the azimuthal period

\begin{equation} T_r = \frac{T_\psi}{2}\,, \end{equation}

Because the radial period is an integer fraction of the azimuthal period, all of the orbits in the homogeneous sphere close and are periodic. From the solved orbit in Equations \(\eqref{eq-homogeneous-orbit-1}\) and \(\eqref{eq-homogeneous-orbit-2}\) it is clear that the azimuthal period is simply \(\omega\) and therefore

\begin{equation} T_\psi = 2\,T_r = \frac{2\pi}{\omega}\,. \end{equation}

The period of every orbit in the homogeneous sphere is therefore the same, the period does not depend on either the energy or angular momentum of the orbit. As we will see below, this is not generally the case. From Equation \(\eqref{eq-tr-tpsi-dpsi}\), we find that \(\Delta \psi\)

\begin{equation}\label{eq-dpsi-homogeneous} \Delta \psi = \pi\,. \end{equation}

Let us look at an example of an orbit in the homogeneous sphere using galpy. As before, we use an isochrone potential with \(b \gg 1\) to give an approximate homogeneous sphere potential at \(r \approx 1\), because galpy does not contain a real homogeneous-sphere potential. We normalize the potential again such that \(v_c(r=1) = 1\). Using this normalization, \(\omega= 1\) and the period \(2\pi/\omega\) in \(x(t)\) and \(y(t)\) should be equal to \(2\pi\). Let’s integrate an orbit with initial condition \((x,y,v_x,v_y) = (1.0,0.0,0.1,1.1)\) in the \(z=0\) plane (note that galpy requires initial conditions to be given in polar coordinates in this case and in general requires initial conditions in cylindrical coordinates). We plot the \(x(t)\) and \(y(t)\) behavior:

In [8]:
from galpy import potential
from galpy.orbit import Orbit
from galpy.util import bovy_plot
hp= potential.IsochronePotential(normalize=1.,b=1000.)
oh= Orbit([1.,0.1,1.1,0.])
ts= numpy.linspace(0.,2.*numpy.pi,1001)
legend(fontsize=18.,loc='upper right',frameon=False);

This demonstrates the sinusoidal behavior expected from Equations \(\eqref{eq-homogeneous-orbit-1}\) and \(\eqref{eq-homogeneous-orbit-2}\). The \(R(t)\) curve should have a period of half of that in \(x(t)\) and \(y(t)\); in this case it should be \(\pi\). This is indeed what we see when we look at the orbit in \(R(t)\):

In [9]:

The orbit in the \((x,y)\) plane is indeed an ellipse with the center at the origin of the coordinate system:

In [10]:
../_images/notebooks_02.-Spherical-Mass-Distributions_26_0.png Orbits in the Kepler potential

Orbits in the homogeneous sphere are about as simple as orbits can get: bodies travel around ellipses at an angular rate that is independent of their energy and angular momentum (or equivalently, of their semi-major axis and eccentricity). Orbits in the Kepler potential are slightly more complicated.

To solve the equations of motion for orbits in the Kepler potential, we use Equation \(\eqref{eq-eom-spher-angle}\) to replace the time derivative with an azimuthal derivative

\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} = \frac{L}{r^2}\,\frac{\mathrm{d}}{\mathrm{d}\psi}\,, \end{equation}

and then we re-write Equation \([\ref{eq-pot-homogeneous}]\) in terms of \(u=1/r\). For general potentials this gives

\begin{equation} \frac{\mathrm{d}^2 u}{\mathrm{d} \psi^2} + u = \frac{1}{L^2\,u^2}\,\frac{\mathrm{d} \Phi}{\mathrm{d} r}\Bigg|_{r=1/u}\,, \end{equation}

and for the specific case of a Kepler potential, this simplifies to

\begin{equation} \frac{\mathrm{d}^2 u}{\mathrm{d} \psi^2} + u = \frac{G\,M}{L^2}\,. \end{equation}

This is the equation of a forced harmonic oscillator with constant forcing, with solution

\begin{equation} u(\psi) = C\,\cos\left(\psi-\psi_0\right)+\frac{G\,M}{L^2}\,. \end{equation}

In terms of \(r(\psi)\) this is

\begin{equation}\label{eq-eom-kepler-orbit} r(\psi) = \frac{1}{C\,\cos\left(\psi-\psi_0\right)+\frac{G\,M}{L^2}}\,. \end{equation}

This is the Equation of an ellipse with the origin at one focus with semi-major axis \(a\) and eccentricity \(e\) given by

\begin{align} a & = \frac{L^2}{G\,M\,(1-e^2)}\,,\\ e & = \frac{C\,L^2}{G\,M}\,, \end{align}

as long as \(e < 1\) (if \(e \geq 1\), the orbit is unbound and described by a parabola for \(e=1\) and a hyperbola for \(e > 1\)).

See BT08 (or other texts on Keplerian orbits) for a derivation of the next set of results. The peri- and apocenter distances are given by

\begin{align} r_p & = a\,(1-e)\,,\\ r_a & = a\,(1+e)\,. \end{align}

The orbital eccentricity \((r_a-r_p)/(r_a+r_p)\) is equal to \(e\) and the orbital eccentricity is therefore equal to the eccentricity of the Keplerian ellipse. This is the motivation for the general definition of the orbital eccentricity in Equation \(\eqref{eq-ecc-spher}\). The energy is given by

\begin{equation}\label{eq-kepler-energy} E = -\frac{G\,M}{2\,a}\,. \end{equation}

The radial period is equal to the azimuthal period and given by

\begin{align} T_r = T_\psi & = 2\pi\,\sqrt{\frac{a^3}{G\,M}}\label{eq-kepler-third-law}\\ & = \frac{\pi}{\sqrt{2}}\,\frac{G\,M}{\sqrt{-E^3}}\,,\label{eq-tr-tpsi-kepler} \end{align}

which is Kepler’s third law. The radial period therefore only depends on the energy, it does not depend on the angular momentum. From Equation \(\eqref{eq-tr-tpsi-dpsi}\), we find that \(\Delta \psi\)

\begin{equation}\label{eq-dpsi-kepler} \Delta \psi = 2\,\pi\,. \end{equation}

Let’s look at the same orbit as we considered above in the homogeneous sphere potential, but now in a point-mass potential. We normalize the potential again such that \(v_c(r=1) = 1\) to get the same normalization as that we used above for the homogeneous sphere. We first compute \(a = (r_p+r_a)/2\) and get the radial period using Equation \(\eqref{eq-tr-tpsi-kepler}\):

In [5]:
kp= potential.KeplerPotential(normalize=1.)
ok= Orbit([1.,0.1,1.1,0.])
rap= ok.rap(analytic=True,pot=kp)
rperi= ok.rperi(analytic=True,pot=kp)
a= 0.5*(rap+rperi)
Tr= 2.*numpy.pi*a**1.5
print("The radial period is %.2f" % Tr)
The radial period is 9.12

and then we integrate the orbit for \(t = T_r\):

In [12]:
ts= numpy.linspace(0.,Tr,1001)
legend(fontsize=18.,loc='upper right',frameon=False);

We see that the orbit is indeed period with a period equal to \(T_r\) also in \(R(t)\):

In [13]:

The orbit covers a much larger range in \(R\) than the same initial conditions in the homogeneous sphere. Finally, let’s compare the orbits in the \((x,y)\) plane:

In [14]:
        label=r'$\mathrm{Point\ mass}$')
oh.plot(label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
legend(fontsize=18.,loc='upper right',frameon=False);

This clearly demonstrates that both orbits are ellipses, but that the origin is at the center of the ellipse for the homogeneous sphere, while it is at the focus of the ellipse for the point-mass potential. Orbits in the isochrone potential

As we discussed already above, the potentials corresponding to a point mass (Kepler) and to a constant density (homogeneous sphere) bracket the range of plausible galactic potentials, as densities typically remain constant or decrease with increasing \(r\). In Sec. 3.2.4, we introduced the isochrone potential, which smoothly interpolates between these two extremes over a scale set by a scale parameter \(b\): for \(r\ll b\), the isochrone potential is approximately equal to that of the homogeneous sphere, at \(r \gg b\), the isochrone potential approaches a Kepler potential. We also mentioned that all orbits in an isochrone potential are analytic, that is, they can be written in terms of elementary functions (like in Equation \(\eqref{eq-homogeneous-orbit-1}\) and Equation \(\eqref{eq-eom-kepler-orbit}\) for the homogeneous sphere and Kepler potential, although the expressions for the isochrone potential are more complicated). Thus, the isochrone potential forms an excellent potential to gain analytical understanding that applies more generally to galactic potentials.

We will not derive the properties of the isochrone potential in detail here, instead refer to Sec. 3.1(c) of BT08. The radial period for an orbit with energy \(E\) and angular momentum \(L\) in the isochrone potential is

\begin{align} T_r & = \frac{\pi}{\sqrt{2}}\,\frac{G\,M}{\sqrt{-E^3}}\,, \end{align}

which is exactly the same expression as for the radial period for an orbit in a Kepler potential (Equation \([\ref{eq-tr-tpsi-kepler}]\)). In an isochrone potential, the radial period also does not depend on \(L\), only on \(E\). The azimuthal range \(\Delta \psi\) that an orbit covers during one radial period is

\begin{equation} \Delta \psi = \pi\,\mathrm{sign}(L)\,\left(1+\frac{|L|}{\sqrt{L^2+4\,G\,M\,b}}\right)\,. \end{equation}

Orbits that never explore \(r \lesssim b\) have \(L \approx r\times v_c \approx r \times \sqrt{G\,M/r}\) or \(L^2 \approx G\,M\,r \gg G\,M\,b\). In this limit, \(\Delta \psi \rightarrow 2\pi\), as required for the Kepler potential (Equation \([\ref{eq-dpsi-kepler}]\)). Orbits that only explore \(r \lesssim b\) are the opposite limit, \(L^2 \ll G\,M\,b\), and therefore \(\Delta \psi \rightarrow \pi\), as required for the homogeneous sphere (Equation \([\ref{eq-dpsi-homogeneous}]\)). The intermediate region smoothly interpolates between these two extremes and we have that

\begin{equation} \pi \leq |\Delta \psi| \leq 2\pi\,, \end{equation}

for the isochrone potential. This in term implies that

\begin{equation}\label{eq-spher-trtphirange} \frac{1}{2} \leq \frac{T_r}{T_\psi} \leq 1\,. \end{equation}

From this, we can conclude the following about the general properties of orbits in galactic potentials:

  • The radial period of orbits is most strongly a function of the energy \(E\) and only a weak function of the angular momentum.
  • During the course of a radial period \(T_r\), orbits perform between half and a full rotation period around the center in \(\psi\). If we approximate the orbit as a precessing, Keplerian ellipse (which only makes sense for mass distributions that are close to Keplerian), the ellipse precesses against the direction of rotation of the body itself.
  • Orbits in galactic potentials are typically not closed: only in the homogeneous sphere or Kepler potential do all orbits close. Orbits instead form rosettes and eventually fill the entire space between \(r_p\) and \(r_a\) as is allowed by their energy and angular momentum.

As an example of the rosette behavior, we integrate an orbit in an isochrone potential with \(b=1\), again normalized sich that \(v_c(r=1) = 1\) and compare it to the orbits we obtained above for the homogeneous sphere and the point-mass potential. To make the rosette behavior more obvious, we start from the same \((x,y)\), but increase the velocity to create a less circular orbit with clear rosette petals:

In [15]:
ip= potential.IsochronePotential(normalize=1.,b=1.)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,4*Tr,1001)
        label=r'$\mathrm{Point\ mass}$')
oh.plot(label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
oi.plot(label=r'$\mathrm{Isochrone}\ (b=1)$',overplot=True)
legend(fontsize=18.,loc='lower left',frameon=False);

The green rosette orbit in the isochrone potential does not close. To see how the orbit moves along this rosette pattern, we generate an animation of the rosette orbit:

In [8]:
ip= potential.IsochronePotential(normalize=1.,b=1.)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,10*Tr,1001)
oi.animate(); # remove the ; to display the animation

End-of-chapter questions

Using dimensional analysis alone, what is a possible value of a hypothetical gravitational potential?
\(1.7\times 10^2\,\mathrm{m\,s}^{-2}\)
IC 2574 is a low-surface-brightness galaxy whose mass budget near the center is dominated by dark matter. Its rotation curve has been measured and found to be linearly rising, \(v_c \propto r\), within about 6 kpc from its center. What type of density profile does that imply for the dark matter at the center of IC 2574?
\(\rho \propto r^{-2}\)
\(\rho \propto r\)
\(\rho \propto r^{-1}\)
\(\rho \propto\) constant
We computed the escape velocity near the Sun assuming that all of the mass of the Milky Way is spherically distributed and contained within the Sun's orbit at 8 kpc from the center. Derive a direct relation between the escape velocity and the circular velocity under this assumption. What is the correct relation?
\(v_\mathrm{esc} = 2\,v_c\)
\(v_\mathrm{esc} = v_c/\sqrt{2}\)
\(v_\mathrm{esc} = \sqrt{3}\,v_c\)
\(v_\mathrm{esc} = \sqrt{2}\,v_c\)
In a homogeneous sphere, the radial period of a circular orbit at 2 kpc from the center is 130 Myr. What is the radial period of the orbit with 62.5% of the energy and half the angular momentum of the circular orbit?
46 Myr
130 Myr
368 Myr
390 Myr
The Sun orbits the center of the Milky Way with a period of approximately 220 Myr. The Sun is not on an exactly circular orbit, but oscillates radially with a small amplitude. Based on what you've learned in this chapter, what is a plausible value for the period of this radial oscillation (note that the answer is not the actual value of this period)?
250 Myr
90 Myr
50 Myr
140 Myr