Warning

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:

(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

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**

and that the potential is related to the density through the **Poisson
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

The second theorem implies that the gravitational potential outside of a shell of radius \(R\) is

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

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

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

or

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\):

or

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

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

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

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

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

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}}\)

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

or

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

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\)

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\))

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

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

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

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

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

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

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

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*,

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

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\)

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\)

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

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
figsize(10,6)
# 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],
xlabel=r'$r$',ylabel=r'$v_c(r)$')
line_homog= potential.plotRotcurve(hp,Rrange=[0.,10.],\
label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
line_plummer= potential.plotRotcurve(pp,Rrange=[0.,10.],\
label=r'$\mathrm{Plummer}\,(b=1)$',overplot=True)
line_isochrone= potential.plotRotcurve(ip,Rrange=[0.,10.],\
label=r'$\mathrm{Isochrone}\,(b=1)$',overplot=True)
legend(handles=[line_kepler[0],line_homog[0],
line_plummer[0],line_isochrone[0]],
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

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

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

and

The circular velocity is given by

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\)

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

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

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

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

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

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
figsize(10,6)
# 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.],\
label=r'$\mathrm{Hernquist}$',yrange=[0.,2.5],
xlabel=r'$r/a$',ylabel=r'$v_c(r/a)$')
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
label=r'$\mathrm{NFW}$',overplot=True)
legend(handles=[line_hernquist[0],line_nfw[0]],
fontsize=18.,loc='upper right',frameon=False)
bovy_plot.bovy_text(r'$\mathrm{Equal\ mass\ at}\ R=12\,a$',
top_left=True,size=18.)
```

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
figsize(10,6)
# 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.],\
label=r'$\mathrm{Hernquist}$',yrange=[0.,2.5],
xlabel=r'$r/a$',ylabel=r'$v_c(r/a)$')
line_nfw= potential.plotRotcurve(nfp,Rrange=[0.,15.],\
label=r'$\mathrm{NFW}$',overplot=True)
legend(handles=[line_hernquist[0],line_nfw[0]],
fontsize=18.,loc='upper right',frameon=False)
bovy_plot.bovy_text(r'$\mathrm{Equal\ inner\ density\ profile}$',
top_left=True,size=18.)
```

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

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

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

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

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

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

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

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

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\)

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

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

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

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

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.

#### 3.3.2.1. 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

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

The equations of motion are therefore

These have the solution

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\):

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

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

The peri- and apocenter radii are simply

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

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

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\)

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
figsize(10,6)
hp= potential.IsochronePotential(normalize=1.,b=1000.)
oh= Orbit([1.,0.1,1.1,0.])
ts= numpy.linspace(0.,2.*numpy.pi,1001)
oh.integrate(ts,hp)
oh.plotx(ylabel=r'$x(t),y(t)$',yrange=[-1.65,1.65],
label=r'$x(t)$')
oh.ploty(label=r'$y(t)$',overplot=True)
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]:
```

```
oh.plotR(yrange=[0.,2.]);
```

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

```
In [10]:
```

```
oh.plot(xrange=[-1.65,1.65],yrange=[-1.65,1.65]);
```

#### 3.3.2.2. 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

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

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

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

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

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

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

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

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

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\)

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)
ok.integrate(ts,kp)
ok.plotx(ylabel=r'$x(t),y(t)$',yrange=[-2.85,2.85],
label=r'$x(t)$')
ok.ploty(label=r'$y(t)$',overplot=True)
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]:
```

```
ok.plotR(yrange=[0.,2.]);
```

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]:
```

```
ok.plot(xrange=[-1.85,1.85],yrange=[-1.85,1.85],
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.

#### 3.3.2.3. 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

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

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

for the isochrone potential. This in term implies that

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)
oi.integrate(ts,ip)
ok.plot(xrange=[-2.55,2.55],yrange=[-2.55,2.55],
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.integrate(ts,ip)
oi.animate(); # remove the ; to display the animation
```

End-of-chapter questions