4.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.
4.2.1. Orbits in the homogeneous sphere¶
The homogeneous sphere (see Section 2.4.2) has particularly simple orbits. First, we re-write the potential for the homogeneous sphere from Equation (2.32) inside of its limiting radius as \begin{equation} \Phi(r) = \frac{1}{2}\,\omega^2\,r^2\,, \end{equation} where we have dropped the constant term and we have written the amplitude in terms of \(\omega^2 = 4\pi\,G\,\rho_0/3\). We only consider orbits that remain at \(r < R\). The equations for the homogeneous sphere are in fact most easily solved in Cartesian coordinates rather than in polar coordinates as in Equation (4.6), so we write this potential in Cartesian coordinates (assuming that the orbital plane is at \(z=0\)): \(\Phi(x,y) = \omega^2\,\left(x^2+y^2\right)/2\). The equations of motion are therefore \begin{align} \ddot{x} &= -\omega^2\,x\,;\quad \ddot{y} = -\omega^2\,y\,. \end{align} These have the solution \begin{align} x(t) &= X\,\cos\left(\omega t + \psi_x\right)\,;\quad y(t) = Y\,\sin\left(\omega t + \psi_y\right)\,. \end{align}
This solution describes an ellipse, with the center at the origin of the mass distribution (\(r=0\)). Note that these are not Keplerian orbits, even though Keplerian orbits are also ellipses; as we will see below, Keplerian orbits are ellipses with the origin of the mass distribution at the focus of the ellipse, not at the center of the ellipse as for the homogeneous sphere. This solution has four constants that need to be determined from the initial condition of the orbit (remember that we are working in the orbital plane, which is constant, and there are therefore only two positions and velocities that count; the position and velocity perpendicular to this plane are zero at all times). These constants are (i+ii) the semi-major and semi-minor axes of the ellipse, \(a\) and \(b\) respectively, (iii) the orientation of the ellipse within the plane, and (iv) the position along the orbit at \(t=0\). Let’s fix the latter two of these by orienting the ellipse such that the semi-major axis is along \(x\) and the orbit is at \(y=0\) at positive \(x\) at \(t=0\): \begin{align}\label{eq-homogeneous-orbit} x(t) & = a\,\cos \omega t\,;\quad y(t) = b\,\sin \omega t\,. \end{align} It is then straightforward to demonstrate that the energy and angular momentum are \begin{align} E & = \frac{\left(a^2+b^2\right)\,\omega^2}{2}\,;\quad L = a\,b\,\omega\,, \end{align} or in terms of the eccentricity \(\varepsilon\) of the ellipse: \(\varepsilon = \sqrt{1-b^2/a^2}\) \begin{align}\label{eq-homogeneous-EL} E & = a^2\,\left(1-\frac{\varepsilon^2}{2}\right)\,\omega^2\,;\quad L = a^2\,\sqrt{1-\varepsilon^2}\,\omega\,. \end{align}
The peri- and apocenter radii are simply \begin{align} r_p & = a\,;\quad r_a = b\,, \end{align} which directly follows from the properties of ellipses. Using Equations (4.22), 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 (4.12) 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 do 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 \(T_r = T_\psi/2\). 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 (4.20), it is clear that the azimuthal period is simply \(\omega\) and therefore \begin{equation} T_\psi = 2\,T_r = 2\pi/\omega\,. \end{equation} The period of every orbit in the homogeneous sphere is therefore the same, the period does not depend on either the energy or angular momentum of the orbit. As we will see below, this is not generally the case. From Equation (4.15), we find that \(\Delta \psi= \pi\).
Let’s look at an example of an orbit in the homogeneous sphere using galpy. We normalize the potential again such that \(v_c(r=1) = 1\) and use a size \(R\) such that the considered orbit never strays outside of the constant-density part of the potential. Using this normalization, \(\omega= 1\) and the period \(2\pi/\omega\) in \(x(t)\) and \(y(t)\) should be equal to \(2\pi\). We 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 Figure 4.2.
[7]:
from galpy import potential
from galpy.orbit import Orbit
figure(figsize=(8,5.5))
hp= potential.HomogeneousSpherePotential(normalize=1.,R=20.)
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)$',gcf=True)
oh.ploty(label=r'$y(t)$',overplot=True)
oh.plotR(label=r'$R(t)$',gcf=True)
legend(fontsize=18.,loc='lower left',frameon=False);
Figure 4.2: An orbit in a homogeneous sphere: \(x(t)\), \(y(t)\), and \(R(t)\).
This figure demonstrates the sinusoidal behavior expected from Equations (4.20). 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 Figure 4.2. Figure 4.3 shows that the orbit in the \((x,y)\) plane is indeed an ellipse with the center at the origin of the coordinate system.
[8]:
oh.plot(xrange=[-1.65,1.65],yrange=[-1.65,1.65])
gca().set_aspect('equal');
Figure 4.3: An orbit in a homogeneous sphere.
As another example, we consider what the orbits of the four inner solar-system planets would look like if the solar system were a homogeneous density rather than having almost all of the mass concentrated in the Sun (as we will consider in the next section). We create a homogeneous sphere density in galpy with a density such that the mass within a sphere of 1 AU is one solar mass (thus, the orbit of the Earth should be similar as in the true potential):
[9]:
import astropy.units as u
# necessary to transfer the solar system's unit system to the potential
from galpy.util.conversion import get_physical
o= Orbit.from_name('solar system')[:4]
hp= potential.HomogeneousSpherePotential(\
amp=2.*numpy.pi/3.*(u.Msun/(4.*numpy.pi*(u.AU)**3./3.)),R=10.*u.AU,
**get_physical(o))
Then we integrate the orbits of the inner solar-system planets for just over one year and plot the resulting orbits (galpy always works in kpc, hence the funny numbers) in Figure 4.4.
[10]:
ts= numpy.linspace(0.,1.1,251)*u.yr
o.integrate(ts,hp)
figure(figsize=(7,7))
o.plot(d1='x',d2='y',xrange=[-1e-8,1e-8],yrange=[-1e-8,1e-8],gcf=True,
label=[rf'$\mathrm{{{planet}}}$'
for planet in Orbit.from_name('solar system').name[:4]])
legend(fontsize=18.,loc='upper right',frameon=False);
Figure 4.4: Orbits in the solar system if the potential were that of a homogeneous sphere.
We see that the orbits are ellipses that all close after one year. The orbit of the Earth is almost circular, just like in the true point-mass potential, but the orbits of the other planets are much more eccentric than their true orbits. Mercury, the innermost planet, has the most eccentric orbit and the most different orbit from its true orbit. This is because the homogeneous density is the most different for Mercury for the way we have set it up.
The following animation shows the orbits of the innermost solar-system planets in the homogeneous density sphere. It clearly demonstrates that the period of all of the orbits is the same:
[11]:
o.animate(d1='x',d2='y'); # remove the ; to display the animation
4.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 (4.7) to replace the time derivative with an azimuthal derivative \begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} = \frac{L}{r^2}\,\frac{\mathrm{d}}{\mathrm{d}\psi}\,, \end{equation} and then we re-write Equation (4.6) in terms of \(u=1/r\). For general potentials this gives \begin{equation}\label{eq-eom-kepler-general-d2udpsi2} \frac{\mathrm{d}^2 u}{\mathrm{d} \psi^2} + u = \frac{1}{L^2\,u^2}\,\frac{\mathrm{d} \Phi}{\mathrm{d} r}\Bigg|_{r=1/u}\,, \end{equation} and for the specific case of a Kepler potential, this simplifies to \begin{equation} \frac{\mathrm{d}^2 u}{\mathrm{d} \psi^2} + u = \frac{G\,M}{L^2}\,. \end{equation} This is the equation of a forced harmonic oscillator with constant forcing, with solution (see Appendix B.4.2) \begin{equation}\label{eq-eom-kepler-orbit} u(\psi) = {1\over r(\psi)} = C\,\cos\left(\psi-\psi_0\right)+\frac{G\,M}{L^2}\,. \end{equation} The combination \(\psi-\psi_0\) is known as the true anomaly. Because \(r\) only depends on \(\psi\) through the cosine, it is clear that orbit closes after one azimuthal period and that \(T_r = T_\psi\). In fact, this is the equation of an ellipse with the origin at one focus with semi-major axis \(a\) and eccentricity \(e\) given by \begin{align}\label{eq-kepler-semimajoraxis} a & = \frac{L^2}{G\,M\,(1-e^2)}\,;\quad e = \frac{C\,L^2}{G\,M}\,, \end{align} as long as \(e < 1\) (if \(e \geq 1\), the orbit is unbound and described by a parabola for \(e=1\) and a hyperbola for \(e > 1\)). That the orbits in the point-mas potential are closed ellipses is known as Kepler’s first law. In terms of the semi-major axis and eccentricity, Equation (4.28) becomes \begin{equation}\label{eq-eom-kepler-orbit-asellipse} r(\psi) = \frac{a\,(1-e^2)}{e\,\cos\left(\psi-\psi_0\right)+1}\,. \end{equation} The peri- and apocenter distances are then obtained by setting \(\psi = \psi_0\) and \(\psi = \pi+\psi_0\): \begin{align} r_p & = a\,(1-e)\,;\quad r_a = a\,(1+e)\,. \end{align} The orbital eccentricity \((r_a-r_p)/(r_a+r_p)\) is equal to \(e\) and the orbital eccentricity is therefore equal to the eccentricity of the Keplerian ellipse. This is the motivation for the general definition of the orbital eccentricity in Equation (4.12).
To obtain the specific energy, we evaluate Equation (4.10) at the pericenter \(r = a(1-e)\), because there \(E = \Phi_\mathrm{eff}(r;L)\) and we can use Equation (4.29) to simplify the resulting expression \begin{align} E & = {L^2 \over 2\,a^2\,(1-e)^2} - {GM \over a\,(1-e)}= -\frac{G\,M}{2\,a}\,.\label{eq-kepler-energy} \end{align}
To compute the azimuthal period, rather than following the general expression of Equation (4.15), we can use Equation (4.7) and (4.30) to get \begin{align} T_\psi & = \int \mathrm{d} t = 2\,\int_{\psi_0}^{\psi_0+\pi} \mathrm{d} \psi\,{\mathrm{d}t \over \mathrm{d} \psi} = {2 \over L}\,\int_{\psi_0}^{\psi_0+\pi} \mathrm{d} \psi\,{r^2(\psi)}\\ & = \frac{2a^2\,(1-e^2)^2}{L}\,\int_{-1}^{1} \mathrm{d} x\,{1 \over \sqrt{1-x^2}\,\left[e\,x+1\right]^2}\\ & = \frac{2\pi\,a^2\,\sqrt{(1-e^2)}}{L}\,. \end{align} As argued above, the radial period \(T_r\) is equal to this azimuthal period. Using Equation (4.29), we can also write this as \begin{equation}\label{eq-kepler-third-law} T_r = T_\psi = \frac{2\pi\,a^{3/2}}{\sqrt{GM}} = = \frac{\pi}{\sqrt{2}}\,\frac{G\,M}{\sqrt{-E^3}}\,, \end{equation} which is Kepler’s third law. The radial period therefore only depends on the energy, it does not depend on the angular momentum.
An alternative parameterization of a Keplerian orbit uses a parameter called the eccentric anomaly \(\eta\) which is related to the true anomaly \(\psi-\psi_0\) through \begin{equation} \sqrt{1-e}\,\tan\left({\psi-\psi_0 \over 2}\right) = \sqrt{1+e}\,\tan\left( {\eta \over 2}\right)\,. \end{equation} In terms of the eccentric anomaly, we can write the orbit as \begin{align}\label{eq-spherpot-kepler-eccanomaly-r} r & = a\,(1-e\,\cos \eta)\,,\\ t & = {T_r \over 2\pi}\,(\eta - e\,\sin \eta)\,.\label{eq-spherpot-kepler-eccanomaly-t} \end{align} The combination \(2\pi t/T_r\) is known as the mean anomaly and Equation (4.39) is Kepler’s equation. Perhaps somewhat surprisingly, the only time we will use this parameterization in this book is when discussing a simple model for the formation of dark-matter halos in Chapter 17.2!
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 (4.36):
[12]:
from galpy import potential
from galpy.orbit import Orbit
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\). Figure 4.5 shows the \(x(t)\), \(y(t)\), and \(R(t)\) behavior.
[13]:
ts= numpy.linspace(0.,Tr,1001)
ok.integrate(ts,kp)
figure(figsize=(8,5.5))
ok.plotx(ylabel=r'$x(t),y(t)$',yrange=[-2.85,2.85],
label=r'$x(t)$',gcf=True)
ok.ploty(label=r'$y(t)$',overplot=True)
ok.plotr(label=r'$R(t)$',gcf=True)
legend(fontsize=18.,loc='lower left',frameon=False);
Figure 4.5: An orbit in a Kepler potential: \(x(t)\), \(y(t)\), and \(R(t)\).
We see that the orbit is indeed periodic with a period equal to \(T_r\). The orbit covers a much larger range in \(r\) than the same initial conditions in the homogeneous sphere. Finally, we compare the orbits for this initial condition in the homogeneous sphere and Kepler potential in the \((x,y)\) plane in Figure 4.6.
[14]:
figure(figsize=(7,7))
ok.plot(xrange=[-2.2,2.2],yrange=[-2.2,2.2],gcf=True,
label=r'$\mathrm{Point\ mass}$')
oh.plot(label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
gca().set_aspect('equal')
legend(fontsize=18.,loc='upper right',frameon=False);
Figure 4.6: An orbit in a Kepler potential compared to a similar orbit in a homogeneous sphere.
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.
As another example, we consider the orbits of planets in the solar system, now in the more realistic point-mass Kepler potential rather than the homogeneous density that we considered in the previous section. We set up a Kepler potential for the solar system:
[15]:
import astropy.units as u
# necessary to transfer the solar system's unit system to the potential
from galpy.util.conversion import get_physical
o= Orbit.from_name('solar system')
kp= potential.KeplerPotential(amp=1.*u.Msun,**get_physical(o))
Then we integrate the orbits for just over one year and plot the resulting orbits of the inner-solar system planets in the left panel of Figure 4.7.
[16]:
figure(figsize=(12,7))
subplot(1,2,1)
ts= numpy.linspace(0.,1.1,251)*u.yr
o.integrate(ts,kp)
o[:4].plot(d1='x',d2='y',xrange=[-1e-8,1e-8],yrange=[-1e-8,1e-8],gcf=True,
label=[rf'$\mathrm{{{planet}}}$'
for planet in Orbit.from_name('solar system').name[:4]])
legend(fontsize=18.,loc='lower left',frameon=False)
gca().set_aspect('equal')
subplot(1,2,2)
ts= numpy.linspace(0.,160.,251)*u.yr
o.integrate(ts,kp)
o[4:].plot(d1='x',d2='y',xrange=[-3e-7,3e-7],yrange=[-3e-7,3e-7],gcf=True,
label=[rf'$\mathrm{{{planet}}}$'
for planet in Orbit.from_name('solar system').name[4:]])
legend(fontsize=18.,loc='lower left',frameon=False)
gca().set_aspect('equal')
tight_layout(w_pad=-0.5)
Figure 4.7: Orbits in the solar system.
We see that in the correct Kepler potential, the orbits of the planets are very close to circular and that, while the orbits of the planets inner to the Earth close after one year, the orbit of Mars does not. This is because the period increases with distance from the Sun according to Kepler’s third law (Equation 4.36). Similarly, we plot the orbits of the outer solar-system planets, integrating for 160 years, in the right panel of Figure 4.7 (we do not plot all eight on the same figure, because the range of scales is so large). The orbit of Neptune is just about to close, because the period of Neptune is about 164 years.
An animation of the orbits of the outer solar-system planets directly shows the increasing period with distance from the Sun:
[17]:
o[4:].animate(d1='x',d2='y'); # remove the ; to display the animation
4.2.3. Orbits in the isochrone potential and other spherical potentials¶
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 Section 2.4.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. 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 orbits in the isochrone potential in detail here, instead leaving this as an exercise. One important property of such orbits is that the radial period for an orbit with specific energy \(E\) and angular momentum \(L\) in the isochrone potential is \begin{align}\label{eq-spherorb-isochrone-Tr} T_r & = \frac{\pi}{\sqrt{2}}\,\frac{G\,M}{\sqrt{-E^3}}\,, \end{align} which is exactly the same expression as for the radial period for an orbit in a Kepler potential (Equation 4.36). In an isochrone potential, the radial period also does not depend on \(L\), only on \(E\), which is the property that gives the isochrone model its name. The azimuthal range \(\Delta \psi\) that an orbit covers during one radial period is \begin{equation}\label{eq-spherorb-isochrone-dpsi} \Delta \psi = \pi\,\mathrm{sign}(L)\,\left(1+\frac{|L|}{\sqrt{L^2+4\,G\,M\,b}}\right)\,. \end{equation} Orbits that never explore \(r \lesssim b\) have \(L \approx r\times v_c \approx r \times \sqrt{G\,M/r}\) or \(L^2 \approx G\,M\,r \gg G\,M\,b\). In this limit, \(\Delta \psi \rightarrow 2\pi\), as required for the Kepler potential, because there \(T_r = T_\psi\). Orbits that only explore \(r \lesssim b\) are in the opposite limit, \(L^2 \ll G\,M\,b\), and therefore \(\Delta \psi \rightarrow \pi\), as required for the homogeneous sphere. The intermediate region smoothly interpolates between these two extremes and we have that \begin{equation} \pi \leq |\Delta \psi| \leq 2\pi\,, \end{equation} for the isochrone potential. This in term implies that \begin{equation}\label{eq-spher-trtphirange} \frac{1}{2} \leq \frac{T_r}{T_\psi} \leq 1\,. \end{equation}
Because over the relatively small radial ranges that most orbits explore, all spherical potentials can be approximated by an isochrone potential for a well-chosen value of \(b\), 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 specific 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 specific energy and angular momentum.
As an example of the rosette behavior, we integrate an orbit in an isochrone potential with \(b=1\), again normalized such 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. The resulting orbit is shown in Figure 4.8.
[18]:
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.,4*Tr,1001)
oi.integrate(ts,ip)
figure(figsize=(7,7))
ok.plot(xrange=[-2.95,2.95],yrange=[-2.95,2.95],gcf=True,
label=r'$\mathrm{Point\ mass}$')
oh.plot(label=r'$\mathrm{Homogeneous\ sphere}$',overplot=True)
oi.plot(label=r'$\mathrm{Isochrone}\ (b=1)$',overplot=True)
gca().set_aspect('equal')
legend(fontsize=18.,loc='lower left',frameon=True,framealpha=0.7,
facecolor='w',edgecolor='w');
Figure 4.8: Orbits in the Kepler, homogeneous sphere, and isochrone potentials.
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 for ten radial periods:
[19]:
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(staticPlot=True); # remove the ; to display the animation
To demonstrate that orbits of the same energy have the same radial period, even though their differing angular momenta can cause them to have quite different orbits, let’s integrate two quite eccentric orbits in an isochrone potential that have the same specific energy, but different angular momenta. The first orbit starts at its apocenter, the second at its pericenter. We also integrate circular orbits at these apo- and pericenter radii to make it clear when each orbit has gone through a radial period.
[20]:
ip= potential.IsochronePotential(normalize=1.,b=0.4)
o1= Orbit([1.,0.,ip.vcirc(1.)-0.7,0.])
os= Orbit([[1.,0.,ip.vcirc(1.)-0.7,0.],
[0.4,0.,numpy.sqrt(2.*(o1.E(pot=ip)-ip(0.4,0.))),0.],
[1.,0.,ip.vcirc(1.),0.], # also animate circular orbits at apo/pericenter
[0.4,0.,ip.vcirc(0.4),0.]])
Tr= o1.Tr(pot=ip)
ts= numpy.linspace(0.,3*Tr,501)
os.integrate(ts,ip)
os.animate(d1=['t','x'],d2=['r','y'],staticPlot=True); # remove the ; to display the animation
A static version of this plot is shown in Figure 4.9.
[21]:
ip= potential.IsochronePotential(normalize=1.,b=0.4)
o1= Orbit([1.,0.,ip.vcirc(1.)-0.7,0.])
os= Orbit([[1.,0.,ip.vcirc(1.)-0.7,0.],
[0.4,0.,numpy.sqrt(2.*(o1.E(pot=ip)-ip(0.4,0.))),0.],
[1.,0.,ip.vcirc(1.),0.], # also animate circular orbits at apo/pericenter
[0.4,0.,ip.vcirc(0.4),0.]])
Tr= o1.Tr(pot=ip)
ts= numpy.linspace(0.,3*Tr,501)
os.integrate(ts,ip)
figure(figsize=(8,4))
subplot(1,2,1)
os.plot(d1='t',d2='r',gcf=True,ls='-')
subplot(1,2,2)
os.plot(d1='x',d2='y',gcf=True,ls='-')
gca().set_aspect('equal')
tight_layout();
Figure 4.9: The radial periods of orbits in the isochrone potential with the same energy, but different angular momentum, are the same.
We see that whenever the blue orbit reaches its apocenter (the green line/circle), the orange orbit reaches its own pericenter (the red line/circle). Thus, their radial periods are the same. However, even though they both start at \(\psi = 0\), at the end of the three displayed radial periods, their azimuthal angle is clearly different (about \(3\pi/2\) for the orange orbit and a little less for the blue orbit), which is due to their different angular momentum.
To illustrate that this behavior is generic (but not exact), let’s consider the same orbits in the spherical NFW potential that represents the dark-matter halo in galpy’s MWPotential2014 model for the Milky Way’s gravitational field. To keep the orbits similar to those in the isochrone example above, we normalize the NFW model also such that \(v_c(r=1) = 1\).
[22]:
from galpy.potential import MWPotential2014
hp= MWPotential2014[2]
# Note that this normalization alters MWPotential2014, so don't do this
# if you still want access to MWPotential2014 in the same Python session
hp.normalize(1.)
o1= Orbit([1.,0.,hp.vcirc(1.)-0.7,0.])
os= Orbit([[1.,0.,hp.vcirc(1.)-0.7,0.],
[0.4,0.,numpy.sqrt(2.*(o1.E(pot=hp)-hp(0.4,0.))),0.],
[1.,0.,hp.vcirc(1.),0.], # also animate circular orbits at apo/pericenter
[0.4,0.,hp.vcirc(0.4),0.]])
Tr= o1.Tr(pot=hp)
ts= numpy.linspace(0.,3*Tr,501)
os.integrate(ts,hp)
os.animate(width=600,d1=['t','x'],d2=['r','y'],staticPlot=True);
A static version of this plot is shown in Figure 4.10.
[23]:
from galpy.potential import MWPotential2014
hp= MWPotential2014[2]
# Note that this normalization alters MWPotential2014, so don't do this
# if you still want access to MWPotential2014 in the same Python session
hp.normalize(1.)
o1= Orbit([1.,0.,hp.vcirc(1.)-0.7,0.])
os= Orbit([[1.,0.,hp.vcirc(1.)-0.7,0.],
[0.4,0.,numpy.sqrt(2.*(o1.E(pot=hp)-hp(0.4,0.))),0.],
[1.,0.,hp.vcirc(1.),0.], # also animate circular orbits at apo/pericenter
[0.4,0.,hp.vcirc(0.4),0.]])
Tr= o1.Tr(pot=hp)
ts= numpy.linspace(0.,3*Tr,501)
os.integrate(ts,hp)
figure(figsize=(8,4))
subplot(1,2,1)
os.plot(d1='t',d2='r',gcf=True,ls='-')
subplot(1,2,2)
os.plot(d1='x',d2='y',gcf=True,ls='-')
gca().set_aspect('equal')
tight_layout();
Figure 4.10: The radial periods of orbits in an NFW potential with the same energy, but different angular momentum, are approximately the same.
It is clear that the radial period of the two orbits is still approximately equal (there is only a small difference after three radial periods of the blue orbit), illustrating that the angular momentum is a minor contributor to the difference in radial period between orbits.
.
.
.