13.2. Orbits in two-dimensional, non-axisymmetric systems¶
To commence our exploration of orbits in triaxial systems, we first consider orbits in a two-dimensional, non-axisymmetric system and introduce a new type of orbits that is important in triaxial systems. We consider a non-axisymmetric version of the logarithmic potential that we first discussed in Chapter 2.4.5 by deforming the potential in the \(y\) direction as \begin{equation}\label{eq-nonaxi-logpot-planar} \Phi(x,y) = \frac{v_c^2}{2}\,\ln\left(x^2+\frac{y^2}{b^2}+R_c^2\right)\,, \end{equation} where in addition to the non-axisymmetry (for \(b\neq 1\)) we have also introduced a core radius \(R_c\) (see also Chapter 8.1.6). At \(x,y \gg R_c\) and \(b \approx 1\), this potential approaches the logarithmic potential from Chapter 2.4.5 in the \(z=0\) plane. As an example, we consider the potential in Equation (13.2) with \(v_c = 1\), \(b=0.9\), and \(R_c = 0.2\). We start by looking at an orbit starting at \((x,y) = (0.1,0)\) with \(v_x = 0\) and \(v_y\) equal to the circular velocity at the initial \(x\) if the potential were axisymmetric.
[10]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
ts= numpy.linspace(0.,30,601)
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
o= Orbit([0.1,0.,lp.vcirc(0.1,phi=0.),0.])
o.integrate(ts,lp)
o.animate(staticPlot=True); # remove the ; to display the animation
A static version of this animation is shown in the left panel of the following figure:
[10]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
ts= numpy.linspace(0.,30,1001)
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
o= Orbit([0.1,0.,lp.vcirc(0.1,phi=0.),0.])
o.integrate(ts,lp)
figure(figsize=(8,4))
subplot(1,2,1)
o.plot(gcf=True)
gca().set_aspect('equal')
subplot(1,2,2)
o.plot(d1='t',d2='Lz',ylabel=r'$L_z$',gcf=True)
ts= numpy.linspace(0.,300,10001)
o.integrate(ts,lp)
axhline(numpy.mean(o.Lz(ts)),color='k',ls='--')
galpy_plot.text(26.,-0.01,r'$\langle L_z \rangle$',size=18.)
tight_layout();
Figure 13.5: A box orbit and its angular momentum as a function of time.
If you play the animation, you see that this orbit does not loop around the center like orbits in an axisymmetric potential, but instead the orbit goes back and forth between clockwise and counterclockwise rotation. It is also clear that this orbit can approach the center arbitrarily closely, again in marked contrast to orbits with non-zero angular momentum in an axisymmetric potential, which are inhibited from doing so by the angular-momentum barrier.
To understand this type of orbit, we consider the potential at small distances from the center. At \(x,y \ll R_c\), the potential becomes \begin{equation}\label{eq-triaxialorbits-2dnonaxi-approxlogpot} \Phi(x,y) \approx \frac{v_c^2}{2\,R_c^2}\,\left[x^2+\frac{y^2}{b^2}\right]\,. \end{equation} This potential can be separated in \(x\) and \(y\) as \(\Phi(x,y) = \Phi_x(x)+\Phi_y(y) = \frac{v_c^2}{2\,R_c^2}\,x^2 + \frac{v_c^2}{2\,R_c^2\,b^2} y^2\); \(\Phi_x(x)\) and \(\Phi_y(y)\) are both the potential of a harmonic oscillator, albeit with different frequencies \(\omega_x = v_c/R_c\) and \(\omega_y = v_c/[R_c\,b]\) in the \(x\) and \(y\) directions, respectively. The motion in \(x\) and that in \(y\) therefore separates into two independent motions in each of \((x,v_x)\) and \((y,v_y)\). Unless \(\omega_x/\omega_y = b\) is a rational number, orbits do not close in the \((x,y)\) plane. Because each of the \(x\) and \(y\) oscillations covers a region between \((-x_\mathrm{max},x_\mathrm{max})\) and \((-y_\mathrm{max},y_\mathrm{max})\), respectively, the orbit in \((x,y)\) covers the entire box between the corners formed by the four combinations of the signs in \((\pm x_\mathrm{max},\pm y_\mathrm{max})\). As an example of this behavior, we look at the same potential as above, but now start an orbit well within the \(x,y \ll R_c\) regime, at \((x,y) = (0.01,0)\), again with \(v_x =0\) and \(v_y\) equal to the circular velocity at \(x\) in the axisymmetric (\(b=1\)) version of this potential. The orbit in \((x,y)\) is displayed in Figure 13.6.
[11]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
ts= numpy.linspace(0.,30,601)
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
o= Orbit([0.01,0.,lp.vcirc(0.01,phi=0.),0.])
o.integrate(ts,lp)
figure(figsize=(4,4))
o.plot(gcf=True)
gca().set_aspect('equal');
Figure 13.6: A box orbit.
Thus, the orbit indeed fills a rectangular box in the \((x,y)\) plane. Note that in this example, \(b\) is a rational number and the orbit of objects at \(x,y \ll R_c\) therefore almost closes. But the anharmonicity resulting from deviations of the true potential from the approximation in Equation (13.3) causes the orbits to not close exactly and slowly fill the box. If you change \(b\) even slightly away from \(0.9\) in the example above, you will see that it fills the box much more quickly (for example, try \(b = \pi/3.5 \approx 0.8976\)).
As we approach \(x,y \approx R_c\), the potential no longer separates into the two harmonic-oscillator components like above, but orbits qualitatively keep following the same behavior: they oscillate between clockwise and counterclockwise rotation—and thus have \(\langle L_z \rangle \approx 0\)—and can approach the center arbitrarily closely. They no longer fill a rectangular box, but still fill a box-like region near the center that fans out at larger \((x,y)\); because of this, all such orbits are known as box orbits. The orbit in Figure 13.5 above is an example of this. The run of the angular momentum as a function of time in the right panel of Figure 13.5 demonstrates that this orbit has no definite sense of rotation; the average value \(\langle L_z \rangle\) of its angular momentum (dashed line) is zero.
When we consider orbits at \(x,y \gg R_c\), we find that the orbits do have a well-defined sense of rotation and loop around the center similar to an orbit in an axisymmetric potential. For example, look at the orbit which has initial \((x,y) = (0.25,0)\) and velocities initialized with a similar procedure as above.
[12]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
ts= numpy.linspace(0.,30,1001)
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
o= Orbit([0.25,0.,lp.vcirc(0.25,phi=0.),0.])
o.integrate(ts,lp)
o.animate(staticPlot=True); # remove the ; to display the animation
This orbit looks like a regular orbit in an axisymmetric potential and is known as a loop orbit. If we consider the sequence of orbits in Figure 13.7 that is initialized in the manner described above, then we see that the orbits go from loop orbits at \(x_0 \gg R_c\) to box orbits at \(x_0 \ll R_c\), with the switch for this family happening around \(x_0 \approx 0.16\). At that point, the loop orbit becomes highly eccentric and turns into a box orbit when \(x_0\) is slightly decreased (third and fourth panel going clockwise).
[12]:
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
ts= numpy.linspace(0.,30,1001)
figure(figsize=(12,8))
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
xinit= [0.25,0.2,0.16,0.15,0.1,0.05]
for ii,xi in enumerate(xinit):
o= Orbit([xi,0.,lp.vcirc(xi,phi=0.),0.])
o.integrate(ts,lp)
if ii > 2: ii= 8-ii # flip bottom row
subplot(2,3,ii+1)
if ii %3 == 0: tylabel= r'$y$'
else: tylabel= None
if ii < 3: txlabel= None
else: txlabel= r'$x$'
o.plot(gcf=True,
xlabel=txlabel,ylabel=tylabel,
xrange=[-.33,0.33],yrange=[-0.33,0.33]);
if ii%3 > 0:
gca().yaxis.set_major_formatter(NullFormatter())
if ii< 3:
gca().xaxis.set_major_formatter(NullFormatter())
galpy_plot.text(rf'$x_0 = {xi:.2f}$',
top_left=True,size=17.)
tight_layout()
Figure 13.7: A sequence of orbits ranging from loopy loops to boxy boxes.
To investigate the orbital structure of this potential, we can again use surfaces of section. For this two-dimensional potential, we can define a surface of section in an analogous way to how we defined it for the meridional plane in an axisymmetric potential above. We can reduce the four-dimensional phase-space from \((x,y,v_x,v_y)\) to \((x,y,v_x)\) by using the conservation of energy (and breaking the sign degeneracy by only keep the positive \(v_y\) solution); then we look at the surface \(y=0\). To systematically study the types of orbits that exist, it is useful to look at surfaces of section of orbits with a fixed specific energy \(E\). Figure 13.8 displays the surface of section for \(E=-0.87\), the energy of the outermost loop orbit above.
[14]:
from scipy import optimize
from galpy.potential import evaluatePotentials
def orbit_xvxE(x,vx,E,pot=None):
"""Returns Orbit at (x,vx,y=0) with given E"""
return Orbit([x,vx,numpy.sqrt(2.*(E-evaluatePotentials(pot,x,0.,phi=0.)
-vx**2./2)),0.])
def zvc(x,E,pot=None):
"""Returns the maximum v_x at this x and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,x,0.,phi=0.)))
lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2)
E= -0.87
x= 0.204
maxvx= zvc(x,E,pot=lp)
vxs= numpy.linspace(0.,maxvx-1e-10,11)
figure(figsize=(8,5),dpi=300)
for vx in vxs:
for sgn in [-1,1]:
o= orbit_xvxE(sgn*x,vx,E,pot=lp)
if vx == vxs[0]: marker= 'o'
else: marker= '.'
o.plotSOS(lp,surface='y',ncross=1_000,
overplot=(vx != vxs[0])+(sgn > 0),
marker=marker,mec='none',zorder=2,
mfc=cm.viridis(numpy.fabs(vx/numpy.amax(vxs))),
xrange=[-0.42,0.42],yrange=[-1.4,1.4],
gcf=True,rasterized=True)
# Also plot the zero-velocity curve, first find min x at which it exists
minx= optimize.brentq(lambda x: E-evaluatePotentials(lp,x,0.,phi=0.),
0.2,0.4)-1e-10
xs= numpy.linspace(-minx,minx,201)
plot(xs,zvc(xs,E,pot=lp),color=cm.viridis(1.),zorder=0)
plot(xs,-zvc(xs,E,pot=lp),color=cm.viridis(1.),zorder=1)
galpy_plot.text(r'$b=0.9$',top_left=True,size=18.);
Figure 13.8: Surface of section of the mildly non-axisymmetric, cored logarithmic potential at \(E = -0.87\).
This surface of section shows the two types of orbits present in the potential. The loop orbits with positive and negative average angular momentum are the loops around \(x\approx \pm 0.204\); the orbit with \((x,v_x) = (0.204,0)\) is the closed loop orbit, which is the parent of this family of loop orbits. The orbits that loop around the center are again box orbits; at constant energy, they circulate around the center in greater loops than the loop orbits. The outermost curve has the maximum amount of energy in the \(x\) component: it has \(v_y \equiv 0\) at all times and thus oscillates back and forth along the \(x\) axis. This orbit is the closed long-axis orbit; it is an extreme version of the box orbits that we have seen in this section and is the parent of the family of box orbits at this energy (analogous to how the closed loop orbit is the parent of the loop orbits). The path of the closed long-axis orbit in the surface of section is a type of zero-velocity curve, because we can obtain the path from the equation \(v_y=0\); we compute and show this as the solid line. We also see that each orbit gives rise to a one-dimensional curve, indicating that each orbit satisfies an integral of the motion beyond the classical energy (a second integral in this case).
Because both the closed loop and the closed long-axis orbit give rise to families of orbits in their vicinity, they are both stable orbits in the sense that a small perturbation to these parent orbits gives a perturbed orbit in the same family as the parent. This is different from the situation for an axisymmetric potential (say, the logarithmic potential that we are considering here with \(b=1\)): the closed long-axis orbit exists in the axisymmetric potential (and has \(L_z=0\)), but any perturbation to it results in a loop orbit, because all orbits with \(L_z \neq 0\) in an axisymmetric potential are loop orbits.
If we increase the level of non-axisymmetry, the importance of box orbits at a given energy increases. Once we go to \(b=0.7\), there are almost no loop orbits left at this energy, as can be seen in Figure 13.9.
[15]:
from scipy import optimize
def orbit_xvxE(x,vx,E,pot=None):
"""Returns Orbit at (x,vx,y=0) with given E"""
return Orbit([x,vx,numpy.sqrt(2.*(E-evaluatePotentials(pot,x,0.,phi=0.)
-vx**2./2)),0.])
def zvc(x,E,pot=None):
"""Returns the maximum v_x at this x and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-evaluatePotentials(pot,x,0.,phi=0.)))
lp= LogarithmicHaloPotential(normalize=True,b=0.7,core=0.2)
E= -0.87
x= 0.0555
maxvx= zvc(x,E,pot=lp)
vxs= numpy.linspace(0.,maxvx-1e-10,11)
figure(figsize=(8,5),dpi=300)
for vx in vxs:
for sgn in [-1,1]:
o= orbit_xvxE(sgn*x,vx,E,pot=lp)
o.plotSOS(lp,surface='y',ncross=1_000,
overplot=(vx != vxs[0])+(sgn > 0),
marker='.',mec='none',zorder=2,
mfc=cm.viridis(numpy.fabs(vx/numpy.amax(vxs))),
xrange=[-0.42,0.42],yrange=[-1.4,1.4],
gcf=True,rasterized=True)
# Also plot the zero-velocity curve, first find min x at which it exists
minx= optimize.brentq(lambda x: E-evaluatePotentials(lp,x,0.,phi=0.),
0.2,0.4)-1e-10
xs= numpy.linspace(-minx,minx,201)
plot(xs,zvc(xs,E,pot=lp),color=cm.viridis(1.),zorder=0)
plot(xs,-zvc(xs,E,pot=lp),color=cm.viridis(1.),zorder=1)
galpy_plot.text(r'$b=0.7$',top_left=True,size=18.);
Figure 13.9: Like Figure 13.8, but for an even more non-axisymmetric potential. Almost all orbits are boxes.
The same behavior occurs if we reduce the energy at constant \(b\): at lower energies, more of the phase-space is taken up by box orbits and the stable closed loop orbits move towards smaller \(|x|\), just like what happens when we decrease \(b\) in Figure 13.9. At the same time, these closed loop orbits become more and more elliptical. At a certain critical energy, the closed loop orbits have \(x=0\) and morph into a closed short-axis orbit. This is expected from the asymptotic small \((x,y)\) behavior that we discussed above: orbits with low energies always remain at \(x,y \ll R_c\) where the potential is approximately that of two independent harmonic oscillators in \(x\) and \(y\); these have closed orbits at \((y,v_y)=(0,0)\), the closed long-axis box orbit, and at \((x,v_x)=(0,0)\), the closed short-axis box orbit. Conversely, when the energy increases, the closed short-axis box orbit turns into two closed loop orbits, one with clockwise and one with counterclockwise rotation. We say that the short-axis box orbit bifurcates into the closed loop orbit.
Loop and box orbits are the main families of orbits in static, two-dimensional non-axisymmetric potentials, with loop orbits the generalization of the stable orbits in axisymmetric potentials and the closed loop orbit the generalization of circular orbits. Before continuing our discussion of orbits in three-dimensional non-axisymmetric potential, we take a detour to learn about another important class of orbits: chaotic orbits.
.
.
.