9.1. Motion in the meridional plane¶
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 3.3 and 3.4, the equations of motion in this coordinate system can be derived by starting from the Lagrangian in cylindrical coordinates \((R,\phi,z)\) given in Equation (3.24) and we then find Equations (3.25). Because \(\phi\) does not explicitly enter the expression for the Lagrangian in Equation (3.24), it is a cyclic coordinate and its conjugate momentum is conserved. 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 forward). 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\). It is then useful to also consider the Hamiltonian from Equation (3.41) with generalized momenta given by Equations (3.39). Because \(p_\phi = \mathrm{constant}\), for an orbit with given \(L_z\), we can write the Hamiltonian 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}= \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 4.1.
The equations of motion in the meridional plane can easily be obtained using Hamilton’s equations applied to the effective Hamiltonian in Equation (9.3): \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 (4.6), 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. Figure 9.1 displays 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).
[6]:
from galpy import potential
from galpy.orbit import Orbit
E, Lz= -1.25, 0.6
ro, vo= 8., 220.
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'
cntrlw= [1.5 for l in levels]
cntrlw[numpy.arange(len(levels))\
[numpy.fabs(levels-E) < 0.01][0]]= 5.
figure(figsize=(10,6))
potential.plotPotentials(potential.MWPotential2014,
rmin=.01,nrs=101,nzs=101,
justcontours=True,
levels=levels,
cntrcolors=cntrcolors,
cntrlw=cntrlw,
effective=True,
Lz=Lz,gcf=True,
ro=ro,vo=vo);
galpy_plot.text(r'$L_z= %.2f$' % Lz,top_left=True,size=18.)
Figure 9.1: The effective potential for an example disk orbit.
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 styled 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 of the effective potential that we introduced for spherical potentials in Chapter 4.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 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 (9.6), 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 (4.6), the solutions of the equations of motion (9.4) and (9.5) 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 (4.12). 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:
[7]:
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
We overplot it on top of the contours of the effective potential for this value of \(L_z\) in the following figure:
[9]:
figure(figsize=(10,6))
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,
cntrlw=cntrlw,
effective=True,
Lz=Lz,gcf=True,
ro=ro,vo=vo)
o.plot(overplot=True,lw=0.5)
galpy_plot.text(5.4,1.9,
r'$E, L_z, v_{R,0}= %.2f, %.2f, %.2f$' % (E,Lz,vR),
size=18.)
Figure 9.2: A disk orbit and its effective potential.
We see that this orbit correctly remains within the region defined by \(E \geq \Phi_\mathrm{eff}(R,z;L_z)\), 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}\). Figure 9.3 instead shows what happens when we start with an orbit that has the same \((E,L_z,R)\), but \(v_R = 0.35\).
[10]:
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)
figure(figsize=(10,6))
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,
cntrlw=cntrlw,
effective=True,
Lz=Lz,gcf=True,
ro=ro,vo=vo)
o.plot(overplot=True,lw=0.5)
galpy_plot.text(5.4,1.9,
r'$E, L_z, v_{R,0}= %.2f, %.2f, %.2f$' % (E,Lz,vR),
size=18.)
Figure 9.3: Another disk orbit and its effective potential.
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. Even integrating the orbits for ten times as long to see whether they ever strays outside of the region defined by their orbits in Figure 9.2 and 9.3, we find that they do not.
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 (9.4) for \(R\) is similar to that for \(r\) in a spherical potential (Equation 4.6), 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 4.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 Figure 9.4.
[11]:
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();
Figure 9.4: Disk orbits displaying rosette patterns.
.
.
.