4.1. General properties of orbits in spherical potentials¶
Equation (4.1), or equivalent forms of it derived in the Lagrangian or Hamiltonian formalism, is an ordinary differential equation that typically needs to be solved numerically. Only for a handful of galactic potentials can Equation (4.1) be solved analytically. Even for most spherical potentials, the full solution of Equation (4.1) needs to be numerical. We will discuss the numerical solution of Equation (4.1) in detail in Section 4.5 below and focus here first on general properties of orbits in spherical potentials specifically.
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}\) as a reminder of this. Because only the radial derivative of the potential is non-zero, the equation of motion in a spherical potential simplifies to \begin{equation}\label{eq-eom-spher} \ddot{\vec{r}} = g_r(r)\,\vec{\hat{r}}\,, \end{equation} where \(\vec{\hat{r}}\) is the radial unit vector and \(\vec{r} = r\,\vec{\hat{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 (Chapter 3.3), this implies the existence of a conserved quantity. This quantity is the (specific) angular momentum vector \begin{equation} \vec{L} = \vec{r}\times \dot{\vec{r}}\,, \end{equation} the cross product of the position vector and the velocity vector. That the angular momentum vector is conserved is most easily demonstrated simply by taking its time derivative. In Chapter 3.1, we showed that the time derivative of the (specific) angular momentum is equal to the (specific) torque \(\vec{N}\) (with a slight abuse of notation, we use the same symbol for torque and specific torque), which is zero for any spherical potential: \begin{align} \dot{\vec{L}} & = \vec{N}= \vec{r}\times\vec{g} = g_r(r)\,r\,\,\vec{\hat{r}}\times\vec{\hat{r}} = 0\,, \end{align} because the cross product of a vector with itself is zero. That the angular momentum vector is conserved implies that the orbit is constrained to the plane through the origin that is perpendicular to the angular momentum vector. We therefore only have to consider the motion in this plane. The magnitude of the angular momentum is \(L\). We use polar coordinates \((r,\psi)\) in this plane to most easily describe the dynamics (using \(\psi\) rather than our regular \(\phi\) as the angle coordinate, as a reminder that this is not the angle \(\phi\) in a cylindrical coordinate frame, rather it is the angle in the orbital plane). We can derive the equations of motion in this coordinate system using the Lagrangian formalism from Chapter 3.3, where the Lagrangian is in this case given by \begin{equation} \mathcal{L} = T-V = \frac{m}{2}\,\left(\dot{r}^2 + r^2\,\dot{\psi}^2\right) - m\Phi(r,\psi)\,. \end{equation} The equations of motion in this plane for \(\Phi(r,\psi) \equiv \Phi(r)\) are found using Lagrange’s equations (3.20) and are given by \begin{align}\label{eq-eom-spher-r} \ddot{r} = \frac{L^2}{r^3} + & g_r(r) = \frac{L^2}{r^3} - \frac{\mathrm{d} \Phi(r)}{\mathrm{d} r}\,\\ \dot{\psi} & = \frac{L}{r^2}\,.\label{eq-eom-spher-angle} \end{align} The second of these equations again simply states the conservation of the angular momentum \(L = r^2\,\dot{\psi}\) and is known as Kepler’s second law in the context of the point-mass potential (but it holds more generally).
The specific energy \(E\) of an orbit in a spherical potential is given by \begin{align} E & = \frac{\dot{r}^2}{2} + \frac{r^2\,\dot{\psi}^2}{2} + \Phi(r)= \frac{\dot{r}^2}{2} + \frac{L^2}{2\,r^2} + \Phi(r)\,\label{eq-energy-spherpot} \end{align} where we have replaced \(\dot{\psi}\) with \(L\) using Equation (4.7). As shown in Section 3.1, this specific energy is conserved if the potential \(\Phi(r)\) does not explicitly depend on time. Because the angular momentum is conserved, the equation of motion for \(r\) in Equation (4.6) can be re-written in terms of an effective potential \(\Phi_\mathrm{eff}(r;L)\) given by \begin{equation} \Phi_\mathrm{eff}(r;L) = \Phi(r) + \frac{L^2}{2\,r^2}\,, \end{equation} where we add \(L\) as an argument to the effective potential to make it explicitly clear that the effective potential is different for orbits with different specific angular momentum. The specific energy from Equation (4.8) is then simply \begin{equation}\label{eq-spherpot-radialenergy} E = \frac{\dot{r}^2}{2} + \Phi_\mathrm{eff}(r;L)\,. \end{equation} The specific 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)\) with \(\Phi(r_\mathrm{max}) = E\), but allows orbits to reach \(r=0\); unbound orbits have \(E > \Phi(\infty)\) and therefore have no \(r_{\mathrm{max}}(E)\) (potentials that do not converge at \(\infty\) therefore only have bound orbits).
The conservation of angular momentum in a spherical potential changes this restriction to \(\Phi_\mathrm{eff}(r;L) \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 at most decreases as fast as the Kepler potential, \(\Phi(r) \propto 1/r\). Therefore, the increase in \(\Phi_\mathrm{eff}(r;L)\) 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 with \(L \neq 0\), the equation \(\Phi_\mathrm{eff}(r;L) = 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 (4.8), we can easily solve for these by solving \(\dot{r}^2 = 0\) \begin{equation}\label{eq-spherpot-drdt} \dot{r}^2 = 2[E-\Phi(r)]-\frac{L^2}{r^2} = 0\,, \end{equation} which is the same as solving \(\Phi_\mathrm{eff}(r;L) = 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 quantify the circularity of an orbit is given by the orbital eccentricity \(e\), which can be defined as \begin{equation}\label{eq-ecc-spher} e = \frac{r_a-r_p}{r_a+r_p}\,. \end{equation} For a circular orbit \(e=0\), while orbits that are close to unbound have \(r_a \gg r_p\) and therefore \(e \rightarrow 1\) as \(r_a\) grows.
The time necessary to travel from pericenter to apocenter and back is known as the radial period \(T_r\). This period can be computed as \begin{align}\label{eq-spherpot-Tr} T_r & = 2\,\int_{r_p}^{r_a}\mathrm{d} t= 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{\mathrm{d} t}{\mathrm{d} r}= 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{1}{\sqrt{2[E-\Phi(r)]-L^2/r^2}}\,, \end{align} where we used Equation (4.11). 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\) (or the rotational period). This is most easily computed by first calculating the change \(\Delta \psi\) in \(\psi\) over a radial period and then computing the number of radial periods necessary to complete one full \(\psi\) rotation. We have that \begin{align}\label{eq-spherpot-dpsi} \Delta \psi & = 2\,\int_{r_p}^{r_a}\mathrm{d} \psi = 2\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{\mathrm{d}t}{\mathrm{d}r}\,\frac{\mathrm{d}\psi}{\mathrm{d}t} = 2\,L\,\int_{r_p}^{r_a}\mathrm{d} r\,\frac{1}{r^2\,\sqrt{2[E-\Phi(r)]-L^2/r^2}}\,, \end{align} where we used Equation (4.11) and Equation (4.7) to substitute \(\mathrm{d}r/\mathrm{d}t\) and \(\mathrm{d}\psi/\mathrm{d}t\), respectively. The azimuthal period is then \begin{equation}\label{eq-tr-tpsi-dpsi} T_\psi = \frac{2\pi}{|\Delta \psi|}\,T_r\,, \end{equation} where we need the absolute value of \(\Delta \psi\), because \(\Delta \psi\) can be negative. All of the quantities \(T_r\), \(\Delta \psi\), and \(T_\psi\) typically need to be computed numerically. Some care is necessary in this numerical calculation, because the integrands integrably diverge at the end-points of the integration interval (see Press et al. 2007 for methods for dealing with this situation). The radial and azimuthal frequencies are simply equal to \(2\pi\) over the radial and azimuthal period.
To convert the description of an orbit in a spherical potential in its orbital plane \((r,\psi)\) to the three-dimensional coordinate system \((r,\phi,\theta)\), we need to consider the orientation of the orbital plane. Because the orbital plane is perpendicular to the angular momentum vector, its orientation is fully specified by the direction of the angular momentum vector. However, conventionally, the orientation of the orbital plane is described by the inclination \(i\), which is the angle between the angular momentum vector and the \(z\) axis, and the \(\phi\) coordinate of the intersection of the orbital plane with the \(z=0\) plane in the quadrant where the orbit passes from negative \(z\) to positive \(z\). This \(\phi\) coordinate is known as the longitude of the ascending node and is somewhat confusingly denoted as \(\Omega\). The inclination is related to the angular momentum vector by \begin{equation}\label{eq-spherorbit-inclination-asLzL} \cos i = {L_z / L}\,. \end{equation} Figure 4.1 demonstrates this geometry graphically.
[6]:
from matplotlib.lines import Line2D
def angle_plot(line1,line2,offset=1,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1):
# Draw angle arc between two lines
# Edited from https://stackoverflow.com/a/25228427 and
# https://gist.github.com/battlecook/0c0bdb7097ec7c8fa160e342b1bf51ef
# Assumes one of the end points == origin
from matplotlib.patches import Arc
# Angle between line1 and x-axis
l1xy= line1.get_xydata()
non_origin= (l1xy[:,0] != origin[0])+(l1xy[:,1] != origin[1])
angle1= numpy.degrees(numpy.arctan2(l1xy[non_origin,1]-origin[1],l1xy[non_origin,0]-origin[0]))[0]
# Angle between line2 and x-axis
l2xy= line2.get_xydata()
non_origin= (l2xy[:,0] != origin[0])+(l2xy[:,1] != origin[1])
angle2= numpy.degrees(numpy.arctan2(l2xy[non_origin,1]-origin[1],l2xy[non_origin,0]-origin[0]))[0]
# Angle between them
theta1= numpy.amin([angle1,angle2])
theta2= numpy.amax([angle1,angle2])
angle= theta2-theta1
return Arc(origin,len_x_axis*offset,len_y_axis*offset,
angle=0.,theta1=theta1,theta2=theta2,
color=color,label=str(angle)+u"\u00b0")
def line_between_points(o1,o2,*args,**kwargs):
return plot([o1[0],o2[0]],[o1[1],o2[1]],*args,**kwargs)
def line_perpendicular(line1,o1,xnew,*args,**kwargs):
l1xy= line1.get_xydata()
slope= -(l1xy[1][0]-l1xy[0][0])/(l1xy[1][1]-l1xy[0][1])
o2= [o1[0]+xnew,o1[1]+slope*xnew]
return line_between_points(o1,o2,*args,**kwargs)
figure(figsize=(6,6))
# We'll draw this as a 2D plot, but with a 3D axes set on it
# so first we draw the axes
xaxis_point= (-0.3,-0.4)
yaxis_point= (-xaxis_point[1]*1.5,0.5*xaxis_point[0])
zaxis_point= (0.,0.75*numpy.sqrt(xaxis_point[0]**2.+xaxis_point[1]**2.))
annotate(text=r'$x$',xy=(0.0,0.0),xytext=xaxis_point,arrowprops=dict(arrowstyle='<|-',color='k'),fontsize=18.,ha='center',va='center')
annotate(text=r'$y$',xy=(0.0,0.0),xytext=yaxis_point,arrowprops=dict(arrowstyle='<|-',color='k'),fontsize=18.,ha='center',va='center')
annotate(text=r'$z$',xy=(0.0,0.0),xytext=zaxis_point,arrowprops=dict(arrowstyle='<|-',color='k'),fontsize=18.,ha='center',va='center')
# dummies for the axes
line_axx= Line2D([0.,xaxis_point[0]],[0.,xaxis_point[1]])
line_axz= Line2D([0.,zaxis_point[0]],[0.,zaxis_point[1]])
# Now we draw the line to the body and the other lines
origin= (0.,0.)
body_vector= [0.37,0.2]
body_vector_xy= [body_vector[0],-0.3]
# Direction of the line of nodes
nodes_vector_xy= numpy.array([-0.06,-0.4])
line_body= line_between_points(origin,body_vector,'k-')[0]
# A dot at the body
plot([body_vector[0]],[body_vector[1]*1.005],'ko',ms=5.)
line_between_points(body_vector,body_vector_xy,'k--')
line_body_xy= line_between_points(origin,body_vector_xy,'k--')[0]
line_nodes= line_between_points(origin,nodes_vector_xy,'k:',lw=2.)[0]
line_orbplane= line_between_points(body_vector,0.6*nodes_vector_xy,'k:',lw=2.)[0]
line_forinc= line_perpendicular(line_nodes,0.6*nodes_vector_xy,0.41,'k-.',lw=2.)[0]
# Indicate angles
gca().add_patch(angle_plot(line_body,line_axz,offset=0.4,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1))
galpy_plot.text(0.13,0.2,r'$\theta$',ha='center',va='center',fontsize=18.)
gca().add_patch(angle_plot(line_axx,line_body_xy,offset=0.575,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1))
galpy_plot.text(0.03,-0.322,r'$\phi$',ha='center',va='center',fontsize=18.)
gca().add_patch(angle_plot(line_nodes,line_body,offset=0.23,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1))
galpy_plot.text(0.15,0.03,r'$\psi$',ha='center',va='center',fontsize=18.)
gca().add_patch(angle_plot(line_axx,line_nodes,offset=0.68,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1))
galpy_plot.text(-0.18,-0.35,r'$\Omega$',ha='center',va='center',fontsize=18.)
gca().add_patch(angle_plot(line_orbplane,line_forinc,offset=0.18,color='k',origin=0.6*nodes_vector_xy,
len_x_axis=1,len_y_axis=1))
galpy_plot.text(0.08,-0.205,r'$i$',ha='center',va='center',fontsize=18.)
# Axis limits etc.
xlim(-0.5,0.5)
ylim(-0.5,0.5)
grid(False)
axis('off');
Figure 4.1: The geometry of an orbit in a spherical potential. A body orbits in the fixed orbital plane defined by the two dotted lines; the angle of the polar coordinates in this plane is \(\psi\). The orbital plane is inclined by an angle \(i\) with respect to the \(x-y\) plane and intersects the \(z=0\) plane in the line-of-nodes (the leftmost dotted line).
We will need this geometry when computing angle coordinates for spherical potentials in Section 4.4 below.
.
.
.