10. Orbits in disks

\label{chapter-diskorbits}

Now that we know how to compute the gravitational potential of a disk-shaped mass distribution, we can start looking at orbits of stars in such potentials. Even though disk galaxies are often characterized by strong non-axisymmetric features like a bar or spiral structure, to a good approximation the mass distributions of these disks are symmetric with respect to rotation around the axis perpendicular to the disk. The other main mass component of galaxies—the dark-matter halo—can also be approximated as begin axisymmetric, at least in the disk plane. When much of the mass in the inner parts of a galaxy is shaped like a bar then that distribution is not well described as an axisymmetric distribution. However, when considering orbits of stars outside of the inner, bar region, the deviation in the gravitational forces between the true bar-shaped mass distribution and an axisymmetric representation are at most a few percent. Therefore, we can still approximate the bar as an axisymmetric distribution (indeed, even as a spherical distribution) for most orbits outside of the bar region. The exception are orbits with specific combinations of radial and azimuthal frequency that make them resonant with the bar’s rotation.

We can therefore usefully investigate disk orbits by considering axisymmetric potentials \(\Phi(R,z;t)\). Moreover, as we discussed in Chapter 2.2.1, the distribution of stars in height above the plane is almost exactly symmetric around \(z=0\), the mid-plane, and therefore \(\Phi(R,z;t) = \Phi(R,-z;t)\). Furthermore, the total mass distribution in the disk region of galaxies changes only very slowly (over many dynamical times) and we can therefore assume that the potential does not depend on time, \(\Phi(R,z;t) \equiv \Phi(R,z)\). We look at the properties of orbits in such potentials in this Chapter.

To illustrate the concepts in this chapter, we will use a simple model for the Milky Way potential, galpy’s MWPotential2014. This is a highly simplified model of the Milky Way’s mass distribution that consists of (i) a Miyamoto-Nagai model (see section 8.2.2) to represent all disk matter (gas plus stars), (ii) an NFW model (see section 3.4.6) to represent the dark-matter halo, and (iii) a spherical bulge potential that is a power-law that is exponentially cut-off at \(r=1.9\,\mathrm{kpc}\) (a model that we did not explicitly discuss, but one can compute its potential and forces using the framework of Section 3.2). This is a convenient model to use because we can compute the total gravitational potential and forces analytically, without doing any of the complicated computations that we discussed in Chapter 8.3.

10.1. Motion in the meridional plane

\label{sec-diskorbits-meridional}

Just as axisymmetric density distributions are best represented in cylindrical coordinates, orbits in these potentials are most easily described in cylindrical coordinates. As already discussed in Chapters 4.3 and 4.4, the equations of motion in this coordinate system can be derived by starting from the Lagrangian in cylindrical coordinates \((R,\phi,z)\)

\begin{equation}\label{eq-lagrangian-cyl} \mathcal{L}(R,\phi,z,\dot{R},\dot{\phi},\dot{z}) = \frac{m}{2}\,\left(\dot{R}^2+[R\,\dot{\phi}]^2+\dot{z}^2\right)-m\,\Phi(R,z)\,. \end{equation}

We can then derive the Hamiltonian by using the momenta \(p_q = \partial \mathcal{L} / \partial \dot{q}\) for \(q=(R,\phi,z)\)

\begin{align} p_R & = \frac{\partial \mathcal{L}}{\partial \dot{R}} = m\,\dot{R}\,,\\ p_\phi & = \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = m\,R^2\,\dot{\phi}\,,\\ p_z & = \frac{\partial \mathcal{L}}{\partial \dot{z}} = m\,\dot{z}\,. \end{align}

Then we get the Hamiltonian

\begin{equation}\label{eq-hamiltonian-cyl} H(R,\phi,z,p_R,p_\phi,p_z) = \frac{1}{2m}\,\left(p_R^2 + \frac{p_\phi^2}{R^2}+p_z^2\right)+m\,\Phi(R,z)\,. \end{equation}

Because \(\phi\) does not explicitly enter the expression for the Lagrangian in Equation \(\eqref{eq-lagrangian-cyl}\), it is a cyclic coordinate and its conjugate momentum is conserved

\begin{equation} \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = m\,R^2\,\dot{\phi} = \mathrm{constant}\,. \end{equation}

We can identify this momentum \(p_\phi\) with the \(z\)-component of the angular momentum \(p_\phi = m\,L_z\), which is therefore conserved (we have introduced the specific angular momentum here and will use this going forwward). Because we are considering time-independent potentials in this section, this means that we have at least two integrals of the motion: the total energy \(E\) and \(L_z\).

Because \(p_\phi = \mathrm{constant}\), for an orbit with given \(L_z\), we can write the Hamiltonian in Equation \(\eqref{eq-hamiltonian-cyl}\) as

\begin{equation} H(R,\phi,z,p_R,p_\phi,p_z) = \frac{1}{2m}\,\left(p_R^2 +p_z^2\right)+m\,\Phi(R,z)+\frac{p_\phi^2}{2m\,R^2}\,. \end{equation}

If we then define the effective potential \(\Phi_\mathrm{eff}(R,z;L_z)\) as

\begin{align} \Phi_\mathrm{eff}(R,z;L_z) & = \Phi(R,z)+\frac{p_\phi^2}{2m^2\,R^2}\nonumber\\ & = \Phi(R,z)+\frac{L_z^2}{2R^2}\,,\label{eq-potential-effective-cyl} \end{align}

orbits with a given \(L_z\) are described by the effective Hamiltonian

\begin{equation}\label{eq-hamiltonian-effective-cyl} H_\mathrm{eff}(R,z,p_R,p_z;L_z) = \frac{1}{2m}\,\left(p_R^2 +p_z^2\right)+m\,\Phi_\mathrm{eff}(R,z;L_z)\,. \end{equation}

Therefore, orbits in a three-dimensional axisymmetric potential are effectively described as a two-dimensional system (with four phase-space dimensions) in \((R,z)\) governed by the effective Hamiltonian. This two-dimensional description is known as the meridional plane. To investigate orbits in axisymmetric potentials, we therefore only have to consider the motion in the meridional plane \((R,z)\) corresponding to the angular momentum \(L_z\) of the orbit. Once we have solved the equations of motion in this plane, we get the time-dependence of \(\phi\) from the conservation of \(L_z\): \(\dot{\phi} = L_z/R^2(t)\). This is analogous to the effective potential that we introduced to describe orbits in spherical potentials in Chapter 5.1.

The equations of motion in the meridional plane can easily be obtained using Hamilton’s equations applied to the effective Hamiltonian in Equation \(\eqref{eq-hamiltonian-effective-cyl}\):

\begin{align}\label{eq-eom-axi-1} {\dot{p_R}\over m} = \ddot{R} & = -\frac{\partial \Phi_\mathrm{eff}(R,z;L_z)}{\partial R} = -\frac{\partial \Phi(R,z)}{\partial R}+\frac{L_z^2}{R^3}\,,\\ {\dot{p_z}\over m} = \ddot{z} & = -\frac{\partial \Phi_\mathrm{eff}(R,z;L_z)}{\partial z} = -\frac{\partial \Phi(R,z)}{\partial z}\,,\label{eq-eom-axi-2} \end{align}

Unfortunately, these equations of motion cannot be solved analytically for essentially any realistic (or not so realistic) disk mass distribution. However, we can solve them numerically. The equation of motion of \(R\) is almost the same as that for the spherical radius in the orbital plane in a spherical potential from Equation (5.7), except that the potential is also a function of \(z\) and not just \(R\).

Before looking at orbits in the meridional plane, let’s take a look at what the effective potential looks like for some disk orbits. Here are contours of the effective potential for \(L_z = 0.6 \times (220\,\mathrm{km\,s}^{-1}) \times (8\,\mathrm{kpc})\) (we express \(L_z\) in units of \((220\,\mathrm{km\,s}^{-1}) \times (8\,\mathrm{kpc})\) such that typical \(L_z\) for stars in the Milky Way disk have \(L_z \approx 1\) in these units):

from galpy import potential
from galpy.orbit import Orbit
E, Lz= -1.25, 0.6
levels= numpy.linspace(-1.25,0.,11)
cntrcolors= ['k' for l in levels]
cntrcolors[numpy.arange(len(levels))\
           [numpy.fabs(levels-E) < 0.01][0]]= '#ff7f0e'
figsize(10,6)
potential.plotPotentials(potential.MWPotential2014,
                         rmin=.01,nrs=101,nzs=101,
                         justcontours=True,
                         levels=levels,
                         cntrcolors=cntrcolors,
                         effective=True,
                         Lz=Lz);
galpy_plot.text(r'$L_z= %.2f$' % Lz,top_left=True,size=18.)

The contours are at \(\Phi_\mathrm{eff} = -1.250\), \(-1.125\), \(-1.00\), \(-0.875\), \(-0.750\), \(-0.625\), \(-0.500\), \(-0.375\), \(-0.250\), \(-0.125\), \(0.00\) (in units of \([220\,\mathrm{km\,s}^{-1}]^2\)) and we have colored the \(\Phi_\mathrm{eff}=-1.250\) contour differently because we will consider orbits with specific energy \(E=-1.250\) below. We have expressed the coordinates in terms of \(R_0 = 8\,\mathrm{kpc}\).

Because energy is conserved and specific energy \(E = (v_R^2+v_z^2)/2+\Phi_\mathrm{eff}(R,z;L_z) \geq \Phi_\mathrm{eff}(R,z;L_z)\), orbits with energy \(E\) are confined to the region \(E \geq \Phi_\mathrm{eff}(R,z;L_z)\). For example, an orbit with \(E=-1.250\) will be confined to the region enclosed by the colored contour above. On this contour, \(E = \Phi_\mathrm{eff}(R,z;L_z)\) and thus \(v_R = v_z = 0\). Therefore, the boundary \(E = \Phi_\mathrm{eff}(R,z;L_z)\) is known as the zero-velocity curve (but note that the rotational velocity \(v_\phi\) is not zero on this curve).

The behavior of the effective potential in the meridional plane is similar to that that we introduced for spherical potentials in Chapter 5.1: As \(R \rightarrow 0\) the term \(L_z^2/(2R^2)\) grows faster than \(\Phi(R,z)\) becomes smaller and the effective potential becomes very large. This is again the angular momentum barrier: an orbit with specific angular momentum \(L_z\) can only reach a non-zero \(R\), no matter how high its energy is. Because we are considering potentials that are symmetric around the \(z=0\) plane, the minimum of the effective potential is at \(z=0\) for all \(R\).

The minimum of the effective potential \(\Phi_\mathrm{eff,min}(L_z)\) is therefore at \((R_\mathrm{min,eff},0)\) for some radius \(R_\mathrm{min,eff}\). The physical meaning of this minimum is clear from the requirement that \(E \geq \Phi_\mathrm{eff}(R,z;L_z)\): an orbit with \(E_c = \Phi_\mathrm{min,eff}(L_z)\) must satisfy \(\Phi_\mathrm{eff}(R,z;L_z) \leq E_c\) or \(\Phi_\mathrm{eff}(R,z;L_z) \leq \Phi_\mathrm{eff,min}(L_z)\) and this orbit can therefore not leave the minimum of the effective potential. Because this minimum occurs at \((R,z) = (R_\mathrm{min,eff},0)\), this orbit therefore has \(R = R_\mathrm{min,eff} = \mathrm{constant}\) and \(z=0\), which is a circular orbit with velocity \(v_c(R_\mathrm{min,eff})\). We can then find \(R_\mathrm{min,eff}\) from the equation of centrifugal balance, noting that \(L_z = R_\mathrm{min,eff}\,v_c(R_\mathrm{min,eff})\):

\begin{equation} \frac{v^2_c(R_\mathrm{min,eff})}{R_\mathrm{min,eff}} = \frac{L_z^2}{R^3_\mathrm{min,eff}}= \frac{\partial \Phi(R,z)}{\partial R}\Bigg|_{R=R_\mathrm{min,eff}}\,.\label{eq-diskorbits-guidingcenter} \end{equation}

This equation needs to be solved numerically, although it can be solved analytically for specific potentials (such as a flattened logarithmic potential, or a power-law potential). It is straightforward to check that one obtains the same equation when solving \(\partial \Phi_\mathrm{eff}(R,z=0;L_z) / \partial R = 0\). The radius of a circular orbit with specific angular momentum \(L_z\) is known as the guiding-center radius and denoted as \(R_g(L_z)\) (the same radius that we denoted as \(R_\mathrm{min,eff}\)). From Equation \eqref{eq-diskorbits-guidingcenter}, it is clear that the guiding-center radius satisfies

\begin{equation} L_z = R_g\,v_c(R_g)\,.\label{eq-diskorbits-guidingcenter-2} \end{equation}

Let’s now take a closer look at orbits in the meridional plane. From the shape of the effective potential and by analogy with the equations of motion in a spherical potential from Equation (5.7), the solutions of the equations of motion \(\eqref{eq-eom-axi-1}\) and \(\eqref{eq-eom-axi-2}\) will be oscillations in \(R\) and \(z\) around the minimum of the effective potential at \((R,z) = (R_\mathrm{min,eff},0)\). Because \(\Phi(R,z)\) for galactic potentials cannot be directly separated into terms involving only \(R\) and only \(z\) (but see below), these oscillations in \(R\) and \(z\) are coupled. Similar to what we did for spherical potentials, we can characterize the orbit with numbers that give a sense of what range in \(R\) and \(z\) the orbit covers: we can define the pericenter \(r_p\) and apocenter \(r_a\) radii as the minimum and maximum of \(\sqrt{R^2+z^2}\) and the orbital eccentricity as \((r_a-r_p)/(r_a+r_p)\), just like for spherical potentials in Equation (5.13). We could also define these based on the cylindrical radius \(R\) rather than based on \(\sqrt{R^2+z^2}\), but for example galpy uses the first definition. To characterize the vertical behavior of the orbit, we define the maximum height above the plane \(z_\mathrm{max}\) that an orbit reaches. It is important to note that these quantities are not as well defined as they were for spherical potentials. As we will see below when we look at example orbits, orbits do not reach \(r_p\), \(r_a\), and \(z_\mathrm{max}\) during every orbit, because of the coupling between the radial and vertical motion, and the numerical determination of these quantities therefore depends on how long you integrate an orbit for. But we will also see that for close-to-circular orbits that stay close to the plane, the radial and vertical motions almost decouple, such that the orbit-to-orbit variations in \(r_p\), \(r_a\), and \(z_\mathrm{max}\) are small for such orbits.

It is useful to look at orbits with the same energy and \(z\) component of the angular momentum, but with different initial \(v_R/v_z\). Therefore, we define a function that returns a galpy orbit with a given \((E,L_z)\) that has initial radial coordinate and velocity \((R,v_R)\) and initial \(z=0\); the initial conditions for the tangential component of the velocity \(v_\phi\) and for the vertical velocity \(v_z\) (chosen to be positive) then follow from the given \((E,L_z)\). We will express these numbers in coordinates where the distance scale is 8 kpc and the velocity scale is \(220\,\mathrm{km\,s}^{-1}\). Let’s start with an orbit with \((E,L_z,R,v_R) = (-1.250,0.6,0.8,0.0)\), integrate it for a dozen dynamical times.

The following animation shows how the orbit evolves in the \((R,z)\) plane:

from galpy import potential
from galpy.orbit import Orbit
def orbit_RvRELz(R,vR,E,Lz,pot=None):
    """Returns Orbit at (R,vR,phi=0,z=0) with given (E,Lz)"""
    return Orbit([R,vR,Lz/R,0.,
                  numpy.sqrt(2.*(E-potential.evaluatePotentials(pot,R,0.)
                                 -(Lz/R)**2./2.-vR**2./2)),0.],ro=8.,vo=220.)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,potential.MWPotential2014)
o.animate(staticPlot=True); # remove the ; to display the animation

Play

Pause

Speed x 2

Speed / 2

Replay

We overplot it on top of the contours of the effective potential for this value of \(L_z\):

potential.plotPotentials(potential.MWPotential2014,
                         rmin=0.3,rmax=1.05,
                         zmin=-0.3,zmax=0.3,
                         nrs=101,nzs=101,
                         justcontours=True,
                         levels=levels,
                         cntrcolors=cntrcolors,
                         effective=True,
                         Lz=Lz)
o.plot(overplot=True,lw=0.5,use_physical=False)
galpy_plot.text(0.675,0.2375,
                r'$E, L_z, v_{R,0}= %.2f, %.2f, %.2f$' % (E,Lz,vR),
                size=18.)

We see that this orbit correctly remains within the region defined by \(E \geq \Phi_\mathrm{eff}(R,z;L_z)\), which in this figure is enclosed by the colored contour, the zero-velocity curve. The orbit touches the zero-velocity curve at four points, but the orbit does not fill the region allowed by \(E \geq \Phi_\mathrm{eff}(R,z;L_z)\). Because we set the initial \(v_R\) to be zero, the initial \(v_z\) is quite high to have energy \(E\) and the orbit reaches quite large heights (the maximum of \(Z\) is \(\approx 0.2 \times 8\,\mathrm{kpc} = 1.6\,\mathrm{kpc}\). Let’s see what happens when we start with an orbit that has the same \((E,L_z,R)\), but \(v_R = 0.35\times 220\,\mathrm{km\,s}^{-1}\):

from galpy import potential
from galpy.orbit import Orbit
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,potential.MWPotential2014)
potential.plotPotentials(potential.MWPotential2014,
                         rmin=0.3,rmax=1.05,
                         zmin=-0.3,zmax=0.3,
                         nrs=101,nzs=101,
                         justcontours=True,
                         levels=levels,
                         cntrcolors=cntrcolors,
                         effective=True,
                         Lz=Lz)
o.plot(overplot=True,lw=0.5,use_physical=False)
galpy_plot.text(0.675,0.2375,
                r'$E, L_z, v_{R,0}= %.2f, %.2f, %.2f$' % (E,Lz,vR),
                size=18.)

Similar to the first orbit, this orbit remains within its zero-velocity curve and touches it at four points, but does not fill the entire region within the zero-velocity curve. Even though these two orbits have the same \(E\) and \(L_z\), they occupy quite different regions of space and therefore must occupy very different regions of phase-space. We can integrate the second orbit for ten times as long to see whether the orbit ever strays outside of the region defined by the orbit above:

from galpy import potential
from galpy.orbit import Orbit
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,1000.,100001)
o.integrate(ts,potential.MWPotential2014)
potential.plotPotentials(potential.MWPotential2014,
                         rmin=0.3,rmax=1.05,
                         zmin=-0.3,zmax=0.3,
                         nrs=101,nzs=101,
                         justcontours=True,
                         levels=levels,
                         cntrcolors=cntrcolors,
                         effective=True,
                         Lz=Lz)
o.plot(overplot=True,ls='None',marker='.',
       ms=.3,alpha=0.3,use_physical=False)
galpy_plot.text(0.675,0.2375,
                r'$E, L_z, v_{R,0}= %.2f, %.2f, %.2f$' % (E,Lz,vR),
                size=18.)

Here we have used transparency to more clearly display where the orbit spends most of its time. We see that even when we integrate the orbit ten times as long as above, it still remains in the same region and does not fill the entire region allowed by its \((E,L_z)\).

So far we have only looked at orbits in the meridional plane \((R,z)\). But while the essential part of their dynamics is contained within this plane, orbits in axisymmetric potentials of course occur in full three-dimensional configuration (or six-dimensional phase space). Therefore, to fully understand orbits in axisymmetric disks, we have to consider the motion in \((x,y)\) in addition to the \((R,z)\) motion. Because the equation of motion \eqref{eq-eom-axi-1} for \(R\) is similar to that for \(r\) in a spherical potential (Equation 5.7), it should come as no surprise that the behavior in the \((x,y)\) plane is similar to that in the orbital plane of a spherical mass distribution. As we saw in Chapter 5.2.3, orbits in \((x,y)\) do not close, but instead form a rosette pattern that eventually fills the entire space between the minimum and maximum radius. The same happens for the projection of a disk orbit in the \((x,y)\) plane, as is clear from the two orbits considered above projected into the \((x,y)\) plane shown in the next figure:

figsize(9,5)
subplot(1,2,1)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,50.,1001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='x',d2='y',gcf=True)
gca().set_aspect(1)
subplot(1,2,2)
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,50.,1001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='x',d2='y',gcf=True)
gca().set_aspect(1)
tight_layout();

10.2. The separability of disk orbits

\label{sec-diskorbits-separable}

The orbits that we looked at above occupy different regions of phase-space, even though they have the same \((E,L_z)\), with the difference being due to their different initial ratio of \(v_R/v_z\). If you play around with changing the \((E,L_z)\) in the examples above, you will see that this is a general property in this Milky-Way-like axisymmetric potential and numerical integrations in a wide range of axisymmetric galactic potentials demonstrate that this is generally the case (Ollongren 1962).

In a completely general, time-dependent potential with no symmetries, orbits would cover a fully six-dimensional region of phase-space if integrated for a long enough time. As we discussed in Chapter 5.3, if a potential has symmetries, Noether’s theorem implies that orbits in this potential satisfy integrals of the motion. For example, as we have already discussed a few times, if the potential does not explicitly depend on time, then the energy of an orbit is conserved (Chapter 4.1). Similarly, when the potential is spherical, the full angular-momentum vector is conserved (as we demonstrated explicitly in Chapter 5.1) and when it is axisymmetric, the component of the angular momentum parallel to the symmetry axis is conserved (see above). As we discussed in Chapter 5.3, one aspect of these integrals that is especially important is that they are all isolating integrals: they restrict the volume of phase-space occupied by an orbit to a lower-dimensional subspace. For example, the conservation of the direction of the angular momentum vector for spherical potential restricts orbits to a plane in configuration space, which reduces the full dimensionality occupied by the orbit in phase-space by two; the conservation of the magnitude of the angular momentum further restrict the motion in the four-dimensional phase-space to a three-dimensional subspace. Conservation of energy for any time-independent potential similarly removes a dimension of phase-space by restricting the orbit to the region where the energy is constant.

If we have identified all integrals of the motion in a given potential, then we expect to see orbits fill the entire region allowed by their integrals of the motion. For example, for general spherical potentials the conservation of energy and the magnitude of the angular momentum imply that orbits are restricted to a radial range \(r_p \leq r \leq r_a\), with the zero-velocity curve in this case given by the two parts: \(r = r_p\) and \(r = r_a\) (“zero velocity” in this case means zero radial velocity; the rotational velocity is always non-zero if the angular momentum is non-zero). Unless the potential is that of a point-mass or a homogeneous sphere (which have an additional isolating integral of motion), the orbit will occupy the entire region \(r_p \leq r \leq r_a\). For example, returning to the example orbit in an isochrone potential with \(b=1\) that we looked at at the end of Chapter 5.2.3 and integrating it for thousands of dynamical times, we see that it eventually covers the entire region between \(r_p \leq r \leq r_a\):

from galpy import potential
from galpy.orbit import Orbit
ip= potential.IsochronePotential(normalize=1.,b=1.)
oi= Orbit([1.,0.25,1.4,0.])
ts= numpy.linspace(0.,10000.,1000001)
oi.integrate(ts,ip)
oi.plot(xrange=[-2.55,2.55],yrange=[-2.55,2.55],
        ls='None',marker='.',
        ms=.3,alpha=0.05,rasterized=True)
gca().set_aspect('equal');

where we have again used transparency to make it clear where the orbit spends most of its time (on the zero velocity curve, similar to what we saw for the axisymmetric orbits above).

The fact that orbits in axisymmetric potentials do not fill the entire region of space allowed by their \((E,L_z)\) indicates that there may be another, third integral beyond \((E,L_z)\) that restricts the motion. It is unclear what this additional integral would be, because the potential does not have any symmetries beyond time independence and axisymmetry. Therefore, we cannot use Noether’s theorem to identify another integral of the motion.

We can gain some insight into what is going on by looking at how well orbits are represented by the product of their motion in the \(z=0\) plane and that perpendicular to this plane. For the axisymmetric orbit with \((E,L_z,v_R) = (0.8,-1.25,0.35)\) that we looked at above, we can plot the projection of the orbit in the \(z=0\) plane and compare it to the orbit that we would get if we integrated the orbit with the same initial \((R,v_R)\) in the \(z=0\) gravitational potential (that is, we completely ignore the vertical motion). We look at the orbit in the \((x,y)\) plane and as \(R\) versus \(t\):

from galpy import potential
from galpy.orbit import Orbit
figsize(10,4.3)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_full= o.plot(d1='x',d2='y',zorder=1,
                  label=r'$\mathrm{Full\ 3D}$',
                  xrange=[-10,12.5],yrange=[-10,12.5],
                  gcf=True)
subplot(1,2,2)
line_full= o.plotR(yrange=[0.,8.],zorder=1,
                  label=r'$\mathrm{Full\ 3D}$',
                  gcf=True)
# Only consider (x,y) part of the orbit
op= o.toPlanar()
op.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_planar= op.plot(d1='x',d2='y',overplot=True,zorder=0,
                    label=r'$\mathrm{Planar\ (2D)}$')
gca().set_aspect('equal')
legend(handles=[line_full[0],line_planar[0]],
       fontsize=18.,loc='upper right',frameon=False)
subplot(1,2,2)
line_planar= op.plotR(overplot=True,zorder=0,
                      label=r'$\mathrm{Planar\ (2D)}$')
legend(handles=[line_full[0],line_planar[0]],
       fontsize=18.,loc='lower right',frameon=False);

We see that we get an excellent approximation to the full orbit’s path in \((x,y)\) by ignoring the vertical part of the orbit. Comparing the time dependence of \(R(t)\) shows that the planar approximation provides a great match to the radial range that is covered by the full orbit. The oscillations in \(R(t)\) are similar in the full orbit and the 2D approximation, but the oscillations start to go out of phase after about a few dozen oscillations. This means that the frequency of the oscillation is not exactly matched in the approximation. But overall we would derive very similar properties about the orbit’s projection on the \(z=0\) plane from approximating it as never leaving the plane. This is remarkable, given that this particular orbit has a rather large eccentricity \(e \approx 0.35\). Most stars, for example, in the Milky-Way disk have smaller eccentricities than this.

We can similarly look at the projection of the orbit on the \((z,v_z)\) plane in phase space and compare it to integrating the same \((z,v_z)\) initial conditions in the fixed one dimensional potential defined by evaluating \(\Phi(R,z)\) at constant \(R\). To get a better sense of how the orbit covers the \((z,v_z)\) plane, we integrate the full orbit for ten times longer than before:

from galpy import potential
from galpy.orbit import Orbit
figsize(10,6)
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.35
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,200.,10001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='z',d2='vz',zorder=0,lw=0.6,
       xrange=[-0.6,0.6],yrange=[-45,45])
# Only consider orbit in z
ol= o.toLinear()
ol.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,R))
line_olmax= ol.plot(overplot=True,zorder=1,
                   label=r'$R = \mathrm{max}\, R(t)$')
olmin= o(ts[numpy.argmin(o.R(ts))]).toLinear()
olmin.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,
             numpy.amin(o.R(ts))))
line_olmin= olmin.plot(overplot=True,zorder=1,
                      label=r'$R = \mathrm{min}\, R(t)$')
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,
             numpy.median(o.R(ts))))
line_olmed= olmed.plot(overplot=True,zorder=1,lw=2.,
                      label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_olmin[0],line_olmed[0],
                line_olmax[0]],
       fontsize=18.,loc='upper right',frameon=False);

The full orbit’s projection is shown in blue. We show three one-dimensional approximations: one initialized at a point along the full orbit that is at the minimum \(R\) that it reaches, one at a point at the median value of \(R\), and one at the maximum \(R\). We see that the 1D approximation to the orbit started at the median \(R\) does a good job of delineating the typical path of the full orbit and the approximations at the minimum and maximum \(R\) capture the extremes of the region well. Nevertheless, the one-dimensional approximation for the \(z(t)\) motion is much less successful than the two-dimensional approximation of the planar part of the motion. This is partly because the orbit has such high eccentricity: the orbit goes through a wide range of \(R\) where the vertical potential is quite different.

This behavior holds more generally. If we look at the first axisymmetric orbit that we considered (the one with \([E,L_z,v_R] = [-1.250,0.6,0.0]\)), we get a similar behavior of the one-dimensional approximations

from galpy import potential
from galpy.orbit import Orbit
R, E, Lz= 0.8, -1.25, 0.6
vR= 0.
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
# 1D little faster if we initialize off of dt=20
ts= numpy.linspace(0.,20.,1001)
o.integrate(ts,potential.MWPotential2014)
ol= o.toLinear()
ol.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,R))
olmin= o(ts[numpy.argmin(o.R(ts))]).toLinear()
olmin.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,
             numpy.amin(o.R(ts))))
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,
             numpy.median(o.R(ts))))
ts= numpy.linspace(0.,200.,10001)
o.integrate(ts,potential.MWPotential2014)
o.plot(d1='z',d2='vz',zorder=0,
       xrange=[-1.74,1.74],yrange=[-165,165.])
line_olmin= olmin.plot(overplot=True,zorder=1,
                      label=r'$R = \mathrm{min}\, R(t)$')
line_olmax= ol.plot(overplot=True,zorder=1,
                   label=r'$R = \mathrm{max}\, R(t)$')
line_olmed= olmed.plot(overplot=True,zorder=1,lw=2.,
                      label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_olmin[0],line_olmed[0],
                line_olmax[0]],
       fontsize=18.,loc='upper right',frameon=False);

It is quite remarkable that we can still delineate the region occupied by the full orbit in the \((z,v_z)\) plane even though this orbit travels up to 1.5 kpc from the plane!

To get a better sense for how this works for more typical disk orbits, which have smaller eccentricities than the orbits that we have looked at so far, we look at the orbit with \((E,L_z,R,v_R) = (-1.385,0.6,0.65,0.02)\). This is orbit is much closer to a circular orbit and has an eccentricity of \(\approx 0.1\). The projection of this orbit onto the \(z=0\) plane and the \(R(t)\) dependence are

from galpy import potential
from galpy.orbit import Orbit
figsize(10,4.3)
R, E, Lz= 0.65, -1.385, 0.6
vR= 0.02
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_full= o.plot(d1='x',d2='y',zorder=1,lw=0.7,
                  label=r'$\mathrm{Full\ 3D}$',
                  xrange=[-7.,9],yrange=[-7.,9.],
                 gcf=True)
subplot(1,2,2)
line_full= o.plotR(yrange=[0.,8.],zorder=1,
                  label=r'$\mathrm{Full\ 3D}$',
                  gcf=True)
op= o.toPlanar()
op.integrate(ts,potential.MWPotential2014)
subplot(1,2,1)
line_planar= op.plot(d1='x',d2='y',overplot=True,zorder=0,
                    lw=0.7,label=r'$\mathrm{Planar\ (2D)}$')
gca().set_aspect('equal')
legend(handles=[line_full[0],line_planar[0]],
       fontsize=18.,loc='upper right',frameon=False)
subplot(1,2,2)
line_planar= op.plotR(overplot=True,zorder=0,
                      label=r'$\mathrm{Planar\ (2D)}$')
legend(handles=[line_full[0],line_planar[0]],
       fontsize=18.,loc='lower right',frameon=False);

where we have again compared the full orbit to the approximation when ignoring the vertical motion. This approximation works similar as before: The radial range covered by the orbit is again well approximated by the 2D approximation, while the oscillation frequency is slightly different, causing the oscillations to drift out of phase as time goes on. Because the orbit covers a smaller range of radii \(R\), the projection of the motion onto the \((z,v_z)\) plane is much better approximated by the one-dimensional approximation at the median \(R\):

from galpy import potential
from galpy.orbit import Orbit
figsize(10,6)
R, E, Lz= 0.65, -1.385, 0.6
vR= 0.02
o= orbit_RvRELz(R,vR,E,Lz,pot=potential.MWPotential2014)
ts= numpy.linspace(0.,20.,1001)
o.integrate(ts,potential.MWPotential2014)
line_full= o.plot(d1='z',d2='vz',zorder=0,lw=0.7,
                 xrange=[-0.3,0.3],yrange=[-35,35.],
                 label=r'$\mathrm{Full\ 3D}$')
ts= numpy.linspace(0.,2.,1001)
olmed= o(ts[o.R(ts) == numpy.median(o.R(ts))]).toLinear()
olmed.integrate(ts,potential.RZToverticalPotential(\
             potential.MWPotential2014,
             numpy.median(o.R(ts))))
line_med= olmed.plot(overplot=True,zorder=1,lw=2.,
                    label=r'$R = \mathrm{med}\, R(t)$')
legend(handles=[line_full[0],line_med[0]],
       fontsize=18.,loc='lower right',frameon=False);

All of these numerical investigations lead us to conclude that a good description of a disk orbit that is close to circular in an axisymmetric potential—a good approximation to the actual dynamical situation of stars in galactic disks—can be obtained by treating the in-plane motion \((R,\phi,v_R,v_\phi)\) and the perpendicular-to-plane motion in \((z,v_z)\) separately. That is, a good approximation of the orbits can be obtained by splitting the overall potential \(\Phi(R,z)\) into a radial \(\Phi_R(R) = \Phi(R,z=0)\) and a vertical \(\Phi_z(z;\tilde{R}) = \Phi(\tilde{R},z)\) part (where the vertical part is defined at some \(\tilde{R}\) along the orbit):

\begin{equation}\label{eq-axipot-sep} \Phi(R,z) = \Phi_R(R) + \Phi_z(z;\tilde{R})\,. \end{equation}

When this is the case, the effective Hamiltonian in Equation \(\eqref{eq-hamiltonian-effective-cyl}\) further simplifies to

\begin{align} H_\mathrm{eff}(R,z,p_R,p_z;L_z) & = \frac{1}{2}\,\left(p_R^2 ++p_z^2\right)+\Phi(R,z)+\frac{L_z^2}{2R^2}\nonumber\\ & = \underbrace{\frac{p_R^2}{2}+\Phi_R(R)+\frac{L_z^2}{2R^2}}_{H_{R,\mathrm{eff}}(R,p_R;L_z)} + \underbrace{\frac{p_z^2}{2}+\Phi_z(z;\tilde{R})}_{H_z(z,p_z)}\label{eq-hamiltonian-effective-sep-cyl}\,. \end{align}

Therefore, in this case, the Hamiltonian splits into two independent parts as \(H_\mathrm{eff}(R,z,p_R,p_z;L_z) = H_{R,\mathrm{eff}}(R,p_R;L_z) + H_z(z,p_z)\), which will therefore have two independent motions.

What this means in terms of integrals of the motion is the following: when the Hamiltonian splits in the manner of Equation \(\eqref{eq-hamiltonian-effective-sep-cyl}\), then the specific energy associated with the radial motion: \(E_R = v_R^2/2 + v_\phi^2/2+\Phi_R(R)\) and that associated with the vertical motion \(E_z = v_z^2 + \Phi(z;\tilde{R})\) are each conserved on their own. Thus, in this case we have three integrals of the motion in this region of phase-space near circular orbits of the axisymmetric potential: \(L_z\), \(E_R\), and \(E_z\). These three integrals of the motion act to restrict the orbit to a three-dimensional subspace of phase-space.

We will use the separability of disk orbits into planar and vertical components to build equilibrium models of galactic disks in Chapter 11. We will come back to the question of whether all orbits in an axisymmetric potential have a third integral in Chapter 14.3, but first we elaborate further on the approximation of the potential in Equation \(\eqref{eq-axipot-sep}\) in more detail.

10.3. Close-to-circular orbits: the epicycle approximation

10.3.1. General considerations

\label{sec-diskorbits-epicycle}

For orbits that are close to circular, we can explicitly derive the separation of the potential in Equation \(\eqref{eq-axipot-sep}\) by expanding the effective potential into a Taylor series around its minimum at the guiding center radius \(R_g\equiv R_g(L_z)\) and \(z=0\):

\begin{equation}\label{eq-diskorbits-epicycleapprox} \Phi_\mathrm{eff}(R,z;L_z) \approx \Phi_\mathrm{eff}(R_g,0;L_z) + \frac{1}{2}\left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial R^2}\right)\Bigg|_{(R_g,0)}\,(R-R_g)^2 +\frac{1}{2}\left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial z^2}\right)\Bigg|_{(R_g,0)}\,z^2\,, \end{equation}

up to terms that are third order in \((R-R_g)\) and \(z\) (the first order terms are zero, because we are at the minimum of the effective potential). The equations of motion \(\eqref{eq-eom-axi-1}\) and \(\eqref{eq-eom-axi-2}\) in this approximation then become

\begin{align} \ddot{R} = \ddot{R}-\ddot{R_g} &= -\left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial R^2}\right)\Bigg|_{(R_g,0)}\,(R-R_g)\,,\\ \ddot{z} & = -\left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial z^2}\right)\Bigg|_{(R_g,0)}\,z\,, \end{align}

(because \(\ddot{R_g} = 0\)). These are the equations of a decoupled harmonic oscillator in \((R-R_g,z)\), with frequencies

\begin{align}\label{eq-eom-epicycle-1} \kappa^2(R_g) & = \left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial R^2}\right)\Bigg|_{(R_g,0)}\,,\\ \nu^2(R_g) & = \left(\frac{\partial^2 \Phi_\mathrm{eff}}{\partial z^2}\right)\Bigg|_{(R_g,0)}\,.\label{eq-eom-epicycle-2} \end{align}

The first of these, \(\kappa\), is known as the epicycle frequency or as the radial frequency, while the second is known as the vertical frequency. A third important frequency is the circular frequency \(\Omega(R_g)\), which is the azimuthal frequency of the circular orbit at \(R_g\) and is therefore \(\Omega(R_g) = v_c(R_g)/R_g = L_z/R_g^2\). Substituting in the equation for the effective potential (Equation \(\ref{eq-potential-effective-cyl}\)), we can express the epicycle and vertical frequencies as

\begin{align} \kappa^2(R_g) & = \left(\frac{\partial^2 \Phi}{\partial R^2}\right)\Bigg|_{(R_g,0)}+3\,\frac{L_z^2}{R_g^4}\,,\\ \nu^2(R_g) & = \left(\frac{\partial^2 \Phi}{\partial z^2}\right)\Bigg|_{(R_g,0)}\,. \end{align}

Using that \(\Omega^2(R_g) = (\partial \Phi/\partial R)|_{R_g}/R_g\), we can also derive the following relation between \(\kappa\) and \(\Omega\)

\begin{equation}\label{eq-epifreqkapp-omega} \kappa^2(R_g) = \left(R\,\frac{\mathrm{d}\Omega^2}{\mathrm{d}R}+4\,\Omega^2\right)\Bigg|_{R_g}\,, \end{equation}

which we can also write as

\begin{equation} \frac{\kappa^2}{\Omega^2}(R_g) = 2\left(\frac{\mathrm{d}\ln [\Omega\,R^2]}{\mathrm{d}\ln R}\right)\Bigg|_{R_g}\,. \end{equation}

As we discussed in Chapter 3.4.4, rotation curves of galaxies go between the extremes of a Keplerian rotation curve \(v_c(R) \propto R^{-1/2}\) and that corresponding to a constant density sphere, \(v_c(R) \propto R\). For these extremes we have that \(\Omega\,R^2 \propto R^{1/2}\) and \(\Omega\,R^2 \propto R^2\), respectively. These correspond to \(\kappa/\Omega = 1\) and \(\kappa/\Omega = 2\) and thus

\begin{equation}\label{eq-epicycle-kappaomegarange} \Omega \lesssim \kappa \lesssim 2\Omega\,, \end{equation}

in galaxies. This is the analogous relation to that for the range in \(T_r/T_\phi\) that we derived for spherical potentials in Equation (5.62). Therefore, we see that the ratio \(\kappa/\Omega\) is informative about the mass distribution.

Because they are the equations of a decoupled harmonic oscillator in \((R-R_g,z)\), it is straightforward to solve the equations of motion \(\eqref{eq-eom-epicycle-1}\) and \(\eqref{eq-eom-epicycle-2}\). The vertical oscillation is

\begin{equation}\label{eq-epicycle-vertical-sol} z = z_\mathrm{max}\,\cos(\nu\,t+\zeta)\,. \end{equation}

where \(z_\mathrm{max}\) is the maximum height above the mid-plane that the orbit reaches and \(\zeta\) is a phase parameter (the phase along the orbit at \(t=0\)). Because disk galaxies have rotation curves that are close to flat, they have the interesting property that the vertical frequency is essentially determined by the local density of matter. To obtain this result, we use the Poisson equation in cylindrical coordinates (Equation 8.25) and substitute \(v_c^2 = R\,\partial \Phi / \partial R\):

\begin{align} 4\pi\,G\,\rho & = \frac{1}{R}\,\frac{\partial}{\partial R}\left(R\,\frac{\partial \Phi}{\partial R}\right) + \frac{\partial^2 \Phi}{\partial z^2}\nonumber\\ & \approx \frac{1}{R}\,\frac{\mathrm{d} v_c^2}{\mathrm{d} R} + \frac{\partial^2 \Phi}{\partial z^2}\nonumber\\ & \approx \nu^2\,, \end{align}

where we use that \(\mathrm{d}v_c / \mathrm{d}R \approx 0\) in galaxies with a flat rotation curve. For the second line to hold at \(z\neq 0\), we need that \(v_c^2 \approx R\,\partial \Phi(R,z) / \partial R\), but this is the case for typical disk mass distributions up to multiple kpc from the mid-plane (shown under very general conditions by Bovy & Tremaine 2012).

Because \(\nu \propto \sqrt{\rho}\), the vertical frequency can only be considered to be a constant when the density is well approximated as being constant. Because galactic disks are thin with a steep vertical density gradient, this clearly cannot be the case for most of the disk. The solution in Equation \(\eqref{eq-epicycle-vertical-sol}\) is therefore only applicable to the small number of objects (e.g., young stars, open clusters) whose orbits do not stray far from the mid-plane (\(z_\mathrm{max} \lesssim 100\,\mathrm{pc}\)).

The motion in \((R-R_g)\) is given by

\begin{equation}\label{eq-epicycle-radial-sol} R(t)-R_g = X\,\cos\left(\kappa\, t + \alpha\right)\,, \end{equation}

with \(X\) and \(\alpha\) constants to be determined from the initial conditions. To obtain the full orbit, we also need to find \(\phi(t)\) from the conservation of angular momentum: \(R^2\,\dot{\phi} = L_z\). Therefore we have that

\begin{equation} \dot{\phi} = \frac{L_z}{R^2(t)} = \frac{L_z}{(R_g+[R-R_g])^2}\,. \end{equation}

Because \(R-R_g \ll R_g\) for a close-to-circular orbit, we can Taylor expand the numerator and get

\begin{equation} \dot{\phi} = \frac{L_z}{R_g^2}\,\left(1-2\,\frac{R(t)-R_g}{R_g}\right)\,. \end{equation}

Substituting in the radial motion from Equation \(\eqref{eq-epicycle-radial-sol}\), that \(\Omega = L_z/R_g^2\), and finally integrating this equation, we find the \(\phi(t)\) solution

\begin{equation}\label{eq-epicycle-angle-sol} \phi(t) = \phi_0+\Omega\,t-2\,\frac{\Omega}{\kappa}\,\frac{X}{R_g}\,\sin\left(\kappa\,t+\alpha\right)\,. \end{equation}

What is the motion that corresponds to the solution in Equations \(\eqref{eq-epicycle-radial-sol}\) and \(\eqref{eq-epicycle-angle-sol}\)? When we subtract out the motion of the guiding center itself, which is \((R,\phi)[t] = (R_g,\Omega\,t+\phi_0)\), we are left with the motion

\begin{align}\label{eq-epicycle-orbit} R(t)-R_g & = X\,\cos\left(\kappa\,t+\alpha\right)\,,\\ R_g\,\left(\phi(t) - \Omega\,t-\phi_0\right) & = -2\,\frac{\Omega}{\kappa}\,X\,\sin\left(\kappa\,t+\alpha\right)\,. \end{align}

This solution corresponds to an ellipse centered on the guiding center with axis ratio given by \(\gamma = 2\Omega/\kappa\). From the range in \(\kappa/\Omega\) in Equation \(\eqref{eq-epicycle-kappaomegarange}\), this axis ratio is \(1/2 \lesssim \gamma \lesssim 1\). This ellipse is known as the epicycle and this explains why we call for \(\kappa\) the epicycle frequency: \(\kappa\) is the frequency at which a star travels around its epicyle ellipse.

10.3.2. The epicycle approximation and the Oort constants

\label{sec-diskorbits-epicycleoort}

We can relate the epicycle approximation to the Oort constants that describe the local velocity field. This allows us to determine how the Oort constants affect the local stellar kinematics beyond the mean field that we discussed in Chapter 9.4.

For an axisymmetric disk, we can relate the frequencies of close-to-circular orbits to the Oort constants. In Equation (9.35), we expressed \(A\) in terms of \(\Omega(R)\). Similarly, from Equation (9.34) for \(B\) and using that \(v_c(R) = \Omega(R)\,R\), we find

\begin{equation}\label{eq-diskorbits-B-asOmega} 2B = -2\Omega(R_0) -R_0\,\frac{\mathrm{d}\Omega}{\mathrm{d}R}\Bigg|_{R_0}\,. \end{equation}

Using Equation \(\eqref{eq-epifreqkapp-omega}\), we can then express \(\kappa(R_0)\) in terms of the Oort constants

\begin{equation}\label{eq-diskorbits-kappa-asOort} \kappa^2(R_0) = -4B\,(A-B)\,. \end{equation}

Using the measured values for \(A\) and \(B\) (Equations 9.37 and 9.38), we therefore find that in the solar neighborhood

\begin{equation} \kappa(R_0) = 33.7\,\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}\,, \end{equation}

and

\begin{equation} \frac{\kappa(R_0)}{\Omega(R_0)} \approx 1.25\,. \end{equation}

Stars in the solar neighborhood therefore execute about 5 radial oscillations for every 4 orbits around the Galactic center.

In Chapter 9.4, we discussed how the mean velocities, as line-of-sight velocities or proper motions, of stars in the solar neighborhood are given in terms of the Oort constants. We can go beyond the mean velocity and determine how the dispersion in velocities is related to the Oort constants when we assume axisymmetry. For that we start from Equation \(\eqref{eq-epicycle-orbit}\) and take the derivative to get expressions for the velocity of a star in the epicyle approximation that is currently located at \(R \approx R_0\)

\begin{align}\label{eq-epicycle-v-1} v_R(t) & = -X\,\kappa\sin\left(\kappa\,t+\alpha\right)\,,\\ v_\phi(t) & = R_0\,\Omega(R_g) -2\,\Omega(R_g)\,X\,\cos\left(\kappa\,t+\alpha\right)\,,\label{eq-epicycle-v-2} \end{align}

where \(X\) is again the amplitude of the star’s epicycle motion in the radial direction, \(R_g\) is its guiding-center radius, and \(R_0\,\Omega(R_g)\) is the contribution to the velocity from the guiding-center’s angular motion. Note that the guiding-center radius \(R_g\) for stars near the Sun is not \(R_0\); in general, each star has its own guiding-center radius in the epicycle approximation.

Assuming that the Sun is on a circular orbit at \(R_0\) (that is, assuming that we have corrected the Sun’s motion to that of the LSR; see Chapter 11.3.1), we can subtract the Sun’s motion from Equations \(\eqref{eq-epicycle-v-1}\) and \(\eqref{eq-epicycle-v-2}\) to get the relative velocity

\begin{align} \Delta v_R(t) = v_R(t) & = -X\,\kappa\sin\left(\kappa\,t+\alpha\right)\,,\\ \Delta v_\phi(t) = v_\phi(t) - v_c(R_0) & = R_0\,\left[\Omega(R_g) - \Omega(R_0)\right]-2\,\Omega(R_g)\,X\,\cos\left(\kappa\,t+\alpha\right)\,, \end{align}

For stars near the Sun, we can Taylor expand \(\Omega(R_g) - \Omega(R_0)\) in terms of \(R_g-R_0\) in the second equation and we only keep terms up to first order in \(X\), which gives

\begin{align} \Delta v_\phi(t) & = R_0\,\frac{\mathrm{d} \Omega}{\mathrm{d} R}\Bigg|_{R_0}\,(R_g-R_0) -2\,\Omega(R_g)\,X\,\cos\left(\kappa\,t+\alpha\right)\nonumber\\ & = -R_0\,\frac{\mathrm{d} \Omega}{\mathrm{d} R}\Bigg|_{R_0}\,X\,\cos\left(\kappa\,t+\alpha\right) -2\,\Omega(R_0)\,X\,\cos\left(\kappa\,t+\alpha\right)\nonumber\\ & = -\left(R_0\,\frac{\mathrm{d} \Omega}{\mathrm{d} R}\Bigg|_{R_0}+2\Omega(R_0)\right)\,X\,\cos\left(\kappa\,t+\alpha\right) \,, \end{align}

where we have used the epicycle approximation for \(R_g-R(t)\) applied to the current location \(R(t) = R_0\). We can then express both components of the relative velocity of local stars measured with respect to the LSR in terms of the Oort constants using Equations \eqref{eq-diskorbits-B-asOmega} and \eqref{eq-diskorbits-kappa-asOort}

\begin{align} \Delta v_R(t) & = -\sqrt{-4B\,(A-B)}\,X\,\sin\left(\kappa\,t+\alpha\right)\,,\\ \Delta v_\phi(t) & = -2\,B\,X\,\cos\left(\kappa\,t+\alpha\right)\,.\label{eq-diskorbits-dvT-epicycle} \end{align}

For a steady-state distribution, we may assume that we are seeing a random mix of phases in a local sample and that the phase and epicycle amplitude are uncorrelated. If the phases are randomly distributed then we have that \(\langle \sin^2\left(\kappa\,t+\alpha\right)\rangle = \langle\cos^2\left(\kappa\,t+\alpha\right)\rangle = 1/2\). Thus, if we average the squared velocities of a sample of local stars, we average over all phases and find the velocity dispersions

\begin{align} \sigma^2_R(t) & = -2\,B\,(A-B)\,\langle X^2\rangle\,,\\ \sigma^2_\phi(t) & = 2\,B^2\,\langle X^2 \rangle\,. \end{align}

If we know the Oort constants, then the radial and tangential velocity dispersion tell us about the typical epicyle amplitude of a sample of stars. For example, intermediate-age stars near the Sun have \(\sigma_R \approx 30\,\mathrm{km\,s}^{-1}\) and therefore have epicycle amplitudes of \(\sqrt{\langle X^2\rangle} \approx \sqrt{\sigma_R^2 / [-2\,B\,(A-B)]}\approx 1.4\,\mathrm{kpc}\).

The ratio of the radial and tangential velocity dispersion of local stars only depends on the Oort constants

\begin{equation}\label{eq-sr-st-epicycle} \frac{\sigma^2_\phi}{\sigma^2_R} = \frac{B}{B-A}\,. \end{equation}

In the solar neighborhood therefore, we expect \(\sigma_\phi/\sigma_R \approx 2/3\) from the measured values of the Oort constants. Because \(B = -A\) for a flat rotation curve, in the case of a flat rotation curve we would find \(\sigma_\phi/\sigma_R = 1/\sqrt{2}\).

Equation \(\eqref{eq-sr-st-epicycle}\) is a useful relation because it shows how the dispersions in the velocities of local stars are directly related to the Galactic potential. However, because for intermediate-age stars the typical epicycle amplitude \(\sqrt{\langle X^2\rangle} \approx 1.4\,\mathrm{kpc}\) is quite large, corrections to Equation \(\eqref{eq-sr-st-epicycle}\) due to non-circular orbits are in fact quite large (Kuijken & Tremaine 1991). Equation \(\eqref{eq-sr-st-epicycle}\) therefore really applies only to the youngest populations of stars near the Sun, but for those stars the assumption of dynamical equilibrium is suspect, because they may not have had enough time to phase mix. To go beyond the description of the velocity dispersion using the epicycle approximation, we need to build equilibrium models for disk populations, which is the subject of the next chapter.