13.3. Chaos in axisymmetric systems¶
The orbit families that we have discussed so far are all regular orbits: they all satisfy three (or two for the planar, non-axisymmetric potentials above) integrals of the motion—including one non-classical one—and they trace simple, smooth one-dimensional curves in the surface of section. Before turning our attention to orbits in triaxial mass distributions, we investigate whether all orbits in three-dimensional, axisymmetric systems or in two-dimensional, non-axisymmetric systems are regular. We will see that they are not: some orbits in some gravitational potentials do not satisfy a non-classical integral of the motion, but rather fill much of the volume of phase-space that is permitted by the conservation of their classical integrals of motion. Such orbits are referred to as ergodic (because for a non-axisymmetric system, only the conservation of energy remains for such orbits and systems for which only the energy matters are ergodic systems; see Chapter 5.6.1) or as chaotic.
To investigate whether all axisymmetric systems satisfy a third integral of the motion, we follow the classical paper by Henon & Heiles (1964), which was the first to demonstrate that this is not the case. As we discussed in Chapter 9.1, motion in an axisymmetric potential can be effectively reduced to the motion in a two-dimensional system in \((R,z)\), the meridional plane, with an effective potential given by the original potential plus a term proportional to the \(z\) component of the angular momentum. Such a system is effectively equivalent to a two-dimensional, non-axisymmetric system if we identify \((R,z) \rightarrow (x,y)\) and continue the potential to negative \(R\) by mirroring the potential \(\Phi(-R,z) = \Phi(R,z)\). Thus, to investigate whether or not an axisymmetric potential has orbits that do not conserve a third integral, we can equivalently study whether or not the effective two-dimensional, non-axisymmetric potential has orbits that do not conserve a second integral.
From this line of argument, Henon & Heiles (1964) further argued that we can study orbits in a simple two-dimensional potential if all we are interested in is an existence proof of ergodic orbits. The potential they proposed is \begin{equation}\label{eq-pot-henon-heiles} \Phi = \frac{1}{2}\,\left(x^2+y^2+2x^2y-\frac{2}{3}y^3\right)= \frac{1}{2}\,\left[R^2+\frac{2\,R^3}{3}\,\sin\left(3\phi\right)\right]\,. \end{equation} One could think of this as the potential in the meridional plane of (a) an orbit with \(L_z = 0\) (such that the \(L_z^2/[2\,R^2]\) term is absent) and (b) for a potential that is not mirror symmetric around the \(z=0\) plane. This is not a realistic model for a galactic potential. However, this simple system exhibits many of the phenomena of chaotic orbits that are important for chaotic orbits in more realistic galactic potential, so it still provides a useful starting point.
We start by considering orbits in the potential in Equation (13.4) with energy \(E=1/12\). Using the techniques of the previous sections, we determine the surface of section for these orbits, this time fixing the \(v_x\) velocity from the conservation of energy and using crossings of the \(x=0\) plane to produce the surface of section in \((y,v_y)\). The surface of section is shown in Figure 13.10.
[53]:
from galpy.potential import HenonHeilesPotential
from galpy import potential
from galpy.orbit import Orbit
from scipy import optimize
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,vx,y,vy) with given E,
with the initial point on the SOS (-sqrt...)"""
return Orbit([y,vy,-numpy.sqrt(2.*(E-potential.evaluateplanarPotentials(pot,y,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-potential.evaluateplanarPotentials(pot,y,phi=numpy.pi/2.)))
figure(figsize=(7,5),dpi=300)
hp= HenonHeilesPotential(amp=1.)
E= 1./12.
# Location in y of closed orbits
ys= [-.30375,0.255,1e-9]
for y in ys:
maxvy= zvc(y,E,pot=hp)
if y > 0.01:
vys= [0.,0.2041/2.,0.2041,maxvy-1e-6] # can't go to maxvy exactly for plotSOS
elif y < -0.01:
vys= [0.,0.2041/2.,0.2041,maxvy-1e-6]
else:
# 0.2041 is location of closed orbit at y=0
vys= [0.2041,-0.2041,0.1641,-0.1641]
for vy in vys:
o= orbit_yvyE(y,vy,E,pot=hp)
if (numpy.fabs(y) > 0.01 and vy == vys[0]) \
or (numpy.fabs(y) < 0.01 and numpy.fabs(numpy.fabs(vy)-0.2041) < 0.01):
marker= 'o'
else: marker= '.'
if y < -0.01: mfc= '#d62728'
elif y > +0.01: mfc= '#ff7f0e'
elif vy > 0: mfc= '#1f77b4'
else: mfc= '#2ca02c'
o.plotSOS(hp,ncross=1_500,surface='x',
overplot=(y != ys[0]) + (vy != vys[0]),
marker=marker,mec='none',zorder=2,mfc=mfc,
xrange=[-0.39,0.53],yrange=[-0.49,0.49],
gcf=True,rasterized=True)
# Add the crossed curve
o= orbit_yvyE(-0.1215,0.,E,pot=hp)
o.plotSOS(hp,ncross=10_000,method='dop853_c',surface='x',
overplot=True,
marker=marker,mec='none',zorder=2,mfc='#9467bd',rasterized=True)
# Also plot the zero-velocity curve
miny= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
-0.4,-0.3)+1e-10
maxy= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
0.4,0.6)-1e-10
ys= numpy.linspace(miny,maxy,201)
plot(ys,zvc(ys,E,pot=hp),color='k',zorder=10,lw=0.5)
plot(ys,-zvc(ys,E,pot=hp),color='k',zorder=11,lw=0.5);
Figure 13.10: Surface of section of the Henon & Heiles (1964) potential at \(E=1/12\).
At \(E=1/12\) there are four stable, closed orbits that are the parents of orbits families; these are indicated by the large dots in Figure 13.10. They span stable orbit families, indicated using the same color as the parents. These four closed orbits are shown in Figure 13.11.
[17]:
figure(figsize=(13,8))
hp= HenonHeilesPotential(amp=1.)
E= 1./12.
# Location in y of closed orbits
ys= [-.30375,0.255,1e-9,1e-9]
vys= [0.,0.,0.2041,-0.2041]
colors= ['#d62728','#ff7f0e','#1f77b4','#2ca02c']
ts= numpy.linspace(0.,30,1001)
for ii,(y,vy) in enumerate(zip(ys,vys)):
o= orbit_yvyE(y,vy,E,pot=hp)
o.integrate(ts,hp,method='dop853_c')
subplot(2,4,ii+1)
if ii == 0: tylabel= r'$y$'
else: tylabel= None
o.plot(gcf=True,
xlabel=None,ylabel=tylabel,color=colors[ii],
xrange=[-.53,0.53],yrange=[-0.53,0.53])
if ii%4 > 0:
gca().yaxis.set_major_formatter(NullFormatter())
gca().xaxis.set_major_formatter(NullFormatter())
if ii < 2:
txtString= rf'$(y_0,v_{{y,0}}) = ({y:.3f},{vy:.1f})$'
else:
txtString= rf'$(y_0,v_{{y,0}}) = ({y:.1f},{vy:.4f})$'
galpy_plot.text(txtString,top_left=True,size=17.)
# Also plot (x,v_y)
subplot(2,4,4+ii+1)
if ii == 0: tylabel= r'$v_y$'
else: tylabel= None
o.plot(d1='x',d2='vy',gcf=True,
xlabel=r'$x$',ylabel=tylabel,color=colors[ii],
xrange=[-.53,0.53],yrange=[-0.43,0.43])
if ii%4 > 0:
gca().yaxis.set_major_formatter(NullFormatter())
tight_layout(w_pad=-0.8)
Figure 13.11: Stable, closed parent orbits of the four main orbit families in the surface of section of Figure 13.10.
The closed orbits with \(v_y=0\) in the surface of section are in the two left columns. These are loop orbits that loop around the origin in a fixed direction. Both of the closed loop orbits follow the same trajectory in \((x,y)\), but they go around this trajectory in opposite directions. This is directly analogous to the circular loop orbits in a two-dimensional axisymmetric potential.
The closed orbits with \(y=0\) in the surface of section are shown in the right two columns of Figure 13.11. They are needle-like, just like the long-axis box orbits above and they are therefore also box orbits.
These four closed orbits are the parents of stable orbit families of loop orbits and of box orbits. For example, if we perturb the initial velocity of these orbits by \(\Delta v_y = 0.05\), we obtain the loop and box orbits displayed in Figure 13.12 that look similar to the loop and box orbits that we discussed in the previous section.
[18]:
figure(figsize=(13,4.5))
hp= HenonHeilesPotential(amp=1.)
E= 1./12.
# Location in y of closed orbits
ys= [-.30375,0.255,1e-9,1e-9]
vys= [0.,0.,0.2041,-0.2041]
colors= ['#d62728','#ff7f0e','#1f77b4','#2ca02c']
# perturb
vys= [vy+0.05 for vy in vys]
ts= numpy.linspace(0.,100,3001)
for ii,(y,vy) in enumerate(zip(ys,vys)):
o= orbit_yvyE(y,vy,E,pot=hp)
o.integrate(ts,hp,method='dop853_c')
subplot(1,4,ii+1)
if ii == 0: tylabel= r'$y$'
else: tylabel= None
o.plot(gcf=True,color=colors[ii],
xlabel=r'$x$',ylabel=tylabel,
xrange=[-.53,0.53],yrange=[-0.53,0.53])
if ii%4 > 0:
gca().yaxis.set_major_formatter(NullFormatter())
gca().xaxis.set_major_formatter(NullFormatter())
if ii < 2:
txtString= rf'$(y_0,v_{{y,0}}) = ({y:.3f},{vy:.2f})$'
else:
txtString= rf'$(y_0,v_{{y,0}}) = ({y:.1f},{vy:.4f})$'
galpy_plot.text(txtString,top_left=True,size=17.)
tight_layout(w_pad=-1.3);
Figure 13.12: Example orbits from the four main orbit families in the surface of section of Figure 13.10.
The purple curve in the surface of section shown in Figure 13.10 straddles the edge of all four stable orbit families. It goes between periods of time where it displays behavior similar to the box orbits above and episodes where it behaves more like the loop orbits above. For example, Figure 13.13 shows four episodes of duration \(\Delta t \approx \mathrm{few\ hundred}\) during a long orbit integration (with total \(\Delta t \approx 13,500\)).
[19]:
tint= [300.,500.,300.,800.]
tinbetween= [3250.,5400.,3000.]
colors= ['#1f77b4','#d62728','#2ca02c','#ff7f0e']
o= orbit_yvyE(-0.1215,0.,E,pot=hp)
figure(figsize=(13,4.5))
for ii in range(4):
ts= numpy.linspace(0.,tint[ii],200001)
o.integrate(ts,hp,method='dop853_c')
subplot(1,4,ii+1)
if ii == 0: tylabel= r'$y$'
else: tylabel= None
o.plot(gcf=True,lw=0.3,
xlabel=r'$x$',ylabel=tylabel,color=colors[ii],
xrange=[-.53,0.53],yrange=[-0.53,0.53])
if ii%4 > 0:
gca().yaxis.set_major_formatter(NullFormatter())
gca().xaxis.set_major_formatter(NullFormatter())
# Go to next loop
if ii < 3:
o= o(ts[-1])
ts= numpy.linspace(0.,tinbetween[ii],200001)
o.integrate(ts,hp,method='dop853_c')
o= o(ts[-1])
tight_layout(w_pad=0.9);
Figure 13.13: Orbit segments along the unstable orbit separating the four main orbit families in the surface of section of Figure 13.10.
In the first, second, and third episodes, the orbit behaves like a member of the box-orbit family, while in the fourth episode, the orbit looks more like a loop orbit. This orbit is unstable: any small perturbation will transform it into a bona-fide member of one of the four stable orbit families. This instability also manifests itself in the numerical orbit integration: even relative energy errors of \(10^{-11}\) cause small changes in the positions and velocities along the orbit that affect how the orbit transitions between the four different behaviors; if we slightly change the parameters of the numerical orbit integration, we get different behavior along the same orbit.
To even better understand this unstable orbit, it is instructive to look at an animation of it. The following animation shows orbit segments over a span of 15,000 time units:
[56]:
o= orbit_yvyE(-0.1215,0.,E,pot=hp)
ts= numpy.linspace(0.,15000.,15001)
o.integrate(ts,hp,method='dop853_c')
o.animate(staticPlot=True); # remove the ; to display the animation
You see that the orbit starts at with the behavior of the positive \(v_y\) box orbit and is stuck there for a long time. It then transitions to the counter-clockwise loop orbit, flips to the clockwise loop orbit, flips back to the counter-clockwise loop orbit, gets caught by the positive \(v_y\) box orbit, goes back to the counter-clockwise loop orbit, goes onto the negative \(v_y\) box orbit, and so on. It’s clear that the behavior is quite erratic, with no clear pattern of when and to what orbit family the orbit will transition at any time.
All the orbits at \(E=1/12\) display one-dimensional curves in the surface of section and they are therefore all regular. However, this system is at the threshold of chaos. If we increase the energy to \(E=1/8\), the surface of section becomes that shown in Figure 13.14.
[54]:
from galpy.potential import HenonHeilesPotential
from galpy import potential
from galpy.orbit import Orbit
from scipy import optimize
def orbit_yvyE(y,vy,E,pot=None):
"""Returns Orbit at (x=0,vx,y,vy) with given E,
with the initial point on the SOS (-sqrt...)"""
return Orbit([y,vy,-numpy.sqrt(2.*(E-potential.evaluateplanarPotentials(pot,y,phi=numpy.pi/2.)
-vy**2./2)),numpy.pi/2.])
def zvc(y,E,pot=None):
"""Returns the maximum v_y at this y and this
energy: the zero-velocity curve"""
return numpy.sqrt(2.*(E-potential.evaluateplanarPotentials(pot,y,phi=numpy.pi/2.)))
figure(figsize=(7,5),dpi=300)
hp= HenonHeilesPotential(amp=1.)
E= 1./8.
# Location in y of closed orbits
ys= [-.37225,0.3025,1e-9,-0.08]
for y in ys:
maxvy= zvc(y,E,pot=hp)
if y > 0.1:
vys= [0.,0.25/8.,0.25/4.,0.113,0.25/2.]
elif y < -0.1:
vys= [0.,0.25/2.,0.2]
elif y < -0.01:
vys= [-0.185,0.1855]
else:
# 0.25 is location of closed orbit at y=0
vys= [0.25,-0.25,0.22,-0.22,-0.19,0.19,0.382]
for vy in vys:
if (numpy.fabs(y) > 0.1 and vy == vys[0]) \
or (numpy.fabs(y) < 0.1 and numpy.fabs(numpy.fabs(vy)-0.25) < 0.01):
marker= 'o'
else: marker= '.'
if y < -0.1: mfc= '#d62728'
elif y > +0.1: mfc= '#ff7f0e'
elif vy > 0.37: mfc= '#d62728' # red island
elif y < -0.01 and vy > 0: mfc= '#1f77b4'
elif y < -0.01 and vy < 0: mfc= '#2ca02c'
elif vy > 0: mfc= '#1f77b4'
else: mfc= '#2ca02c'
o= orbit_yvyE(y,vy,E,pot=hp)
o.plotSOS(hp,ncross=1_500,surface='x',
overplot=(y != ys[0]) + (vy != vys[0]),
marker=marker,mec='none',zorder=2,mfc=mfc,
xrange=[-0.47,0.73],yrange=[-0.53,0.53],
gcf=True,rasterized=True)
# Add the chaotic regions, need to integrate longer
oc= orbit_yvyE(-0.1,0.,E,pot=hp)
oc.plotSOS(hp,ncross=20_000,surface='x',
overplot=True,
marker='.',ms=2.,mec='none',zorder=2,mfc='#9467bd',rasterized=True)
# closed orbit in the chaotic regime
o= orbit_yvyE(-0.12,0.,E,pot=hp)
o.plotSOS(hp,ncross=1000,surface='x',
overplot=True,
marker='o',ms=7.,mec='none',zorder=2,mfc='#e377c2',rasterized=True)
# Also plot the zero-velocity curve
miny= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
-0.45,-0.35)+1e-10
maxy= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
0.6,0.8)-1e-10
ys= numpy.linspace(miny,maxy,201)
plot(ys,zvc(ys,E,pot=hp),color='k',zorder=10,lw=0.5)
plot(ys,-zvc(ys,E,pot=hp),color='k',zorder=11,lw=0.5);
Figure 13.14: Surface of section of the Henon & Heiles (1964) potential at \(E=1/8\).
This surface of section looks very different from the one above. The same four stable orbit families are still present, but they are now embedded in a sea of purple points, all of which are points from a single orbit. This orbit does not follow a simple one-dimensional curve, but in fact fills an two-dimensional area. It therefore does not satisfy a second integral of the motion in addition to the energy and we have found a non-regular orbit (or a chaotic or ergodic orbit).
In addition to the chaotic orbit, Figure 13.14 also displays interesting structure around the stable orbit families. The stable orbit families themselves look similar to the orbit families at the \(E=1/12\) energy at which all orbits are regular: there are four closed parent orbits at the center of these families, indicated by the dots; these closed orbits look similar to those for \(E=1/12\) above. The boundaries of the orbit families no longer touch, but are now separated by the chaotic region. However, we see that there are small areas around the stable orbit families that the chaotic orbits do not enter. These are so-called islands, surrounding the main landmass of the stable orbit family. The loop-orbit families are surrounded by five islands, while the box-orbit families have four. The approximate edge of these islands are indicated with points with the same color as the main orbit family in the figure above. All islands surrounding a given family are visited by a single orbit: an orbit jumps between the different islands of the same color over the course of the orbit integration. The centers of the islands are higher-order closed orbits: they pierce through the surface of section at five (for the loop orbits) and four (for the box orbits) distinct points, located at the centers of the islands.
The islands of the orange loop-orbit family display further curious behavior, because the edge locus does not fully close. This is because we have in fact not found the actual edge of the island, but have instead landed on a higher-order set of islands. In Figure 13.15, we zoom in on two of the orange islands, the top-left and the top-middle one from Figure 13.14.
[55]:
figure(figsize=(12,5),dpi=300)
ys= [0.3025,0.3025,0.3025,0.3025,0.35,0.35,0.3525]
vys= [0.125,0.113,0.13,0.14,0.14,0.1325,0.15]
for y,vy in zip(ys,vys):
o= orbit_yvyE(y,vy,E,pot=hp)
sectys,sectvys= o.SOS(hp,ncross=1_500,surface='x')
subplot(1,2,1)
if vy == vys[-1]: marker=' o'
else: marker= '.'
if numpy.fabs(vy-0.14) < 0.001 \
and numpy.fabs(y-0.3025) < 0.001: mfc= '#2ca02c'
else: mfc= '#ff7f0e'
plot(sectys,sectvys,marker,mec='none',zorder=2,mfc=mfc,ms=5.,rasterized=True)
subplot(1,2,2)
plot(sectys,sectvys,marker,mec='none',zorder=2,mfc=mfc,ms=5.,rasterized=True)
# Overplot chaotic trajectory from above
subplot(1,2,1)
xmin, xmax= -0.06,0.16
ymin, ymax= 0.03,0.12
sectys,sectvys= oc.SOS(hp,ncross=20_000,surface='x')
indx= (sectys > xmin)*(sectys < xmax)\
*(sectvys > ymin)*(sectvys < ymax)
plot(sectys[indx],sectvys[indx],'.',mec='none',zorder=2,mfc='#9467bd',rasterized=True)
xlabel(r'$y$')
ylabel(r'$v_y$')
xlim(xmin,xmax)
ylim(ymin,ymax)
subplot(1,2,2)
xmin, xmax= 0.26,0.44
ymin, ymax= 0.1,0.19
indx= (sectys > xmin)*(sectys < xmax)\
*(sectvys > ymin)*(sectvys < ymax)
plot(sectys[indx],sectvys[indx],'.',mec='none',zorder=2,mfc='#9467bd',rasterized=True)
xlabel(r'$y$')
xlim(xmin,xmax)
ylim(ymin,ymax)
tight_layout();
Figure 13.15: Close-up of the surface of section of Figure 13.14 near some of the islands.
The lower orange edge is the edge of the main loop-orbit family and the orange patches are the loci that are displayed in the full surface of section; we have included some more curves closer to the center for each island as well. We see that in this zoomed-in version the loci near the edge of each island in the full surface of section turn out to be a family of eight islands themselves. These islands are second-order islands: they are islands of islands. At the center of these islands is another, more complicated, closed orbit. If we were to zoom in on these islands, we would find further hierarchies of islands. We have also found another set of islands inside of the main island; these are indicated in green and illustrate how ubiquitous islands are in the surface of section. The deep hierarchy of islands of ever decreasing size around islands that are themselves islands, etc., leads the chaos in the following way: as we go to orbits that are part of a set of islands that is very deep down in the hierarchy of islands, the orbit will appear to jump randomly between different parts of the surface of section and will appear to densely fill a two-dimensional area.
There is one more orbit shown in the full surface of section in Figure 13.14 that is worth commenting on: the orbit that gives rise to the pink points. This orbit is a periodic orbit that appears to be deep in the chaotic regime. Any small perturbation to this orbit turns it into a chaotic orbit. For example, Figure 13.16 shows this period orbit over \(\Delta t = 10^5\) together with an orbit with a slightly lower initial \(y\).
[60]:
figure(figsize=(10,5),dpi=250)
subplot(1,2,1)
ts= numpy.linspace(0.,100000.,3000001)
o= orbit_yvyE(-0.12,0.,E,pot=hp)
o.integrate(ts,hp)
o.plot(gcf=True,lw=0.5,xrange=[-0.7,0.7],yrange=[-0.55,0.75],rasterized=True)
galpy_plot.text(r'$y_0=-0.120$',top_left=True,size=18.)
subplot(1,2,2)
# Integrate for same length, but get denser sampling
ts= numpy.linspace(0.,100000.,10000001)
o= orbit_yvyE(-0.125,0.,E,pot=hp)
o.integrate(ts,hp)
o.plot(gcf=True,marker='.',ms=0.3,ls='none',alpha=0.013,rasterized=True,
xrange=[-0.7,0.7],yrange=[-0.55,0.75])
galpy_plot.text(r'$y_0=-0.125$',top_left=True,size=18.)
tight_layout();
Figure 13.16: Chaotic orbit (right) and a closed orbit near it (left) in the surface of section of Figure 13.14.
These orbits look very different: the periodic orbit remains well-behaved for long times, while the chaotic orbit fills all of the phase-space allowed by the conservation of energy (we have used transparency to more clearly display the projection of this orbit in \([x,y]\)).
The discussion in this section makes it clear that not all orbits in every two-dimensional non-axisymmetric potential have a second integral and consequently, that not all orbits in every three-dimensional axisymmetric potential satisfy a third integral. Some realistic galactic potentials consisting of a disk and halo also have this property (e.g., Hunter 2005), but not all do. For some axisymmetric potentials all orbits are regular. An example of this is the Kuzmin potential that we discussed in Chapter 7.2.1. Such potentials are said to be integrable. Whether or not a galactic potential contains chaotic orbits and what fraction of the phase-space is chaotic is a question that needs to be answered for each potential separately.
It should be clear from this section how useful surfaces of section are in studying the importance of chaos in gravitational potentials. Because chaotic orbits fill a two-dimensional area in the surface of section, but do not enter regions of stability, the surface of section of a single chaotic orbit immediately delineates the chaotic and possibly-non-chaotic regions of phase space. For example, if we just plot the surface of section of the chaotic orbit at \(E=1/8\) from above, we get the result in Figure 13.17.
[52]:
figure(figsize=(7,5),dpi=300)
E= 1./8.
oc.plotSOS(hp,ncross=40_000,surface='x',gcf=True,
marker='.',ms=2.,mec='none',zorder=2,mfc='#9467bd',rasterized=True)
# Also plot the zero-velocity curve
miny= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
-0.45,-0.35)+1e-10
maxy= optimize.brentq(lambda y: E-potential.evaluateplanarPotentials(hp,y,phi=numpy.pi/2.),
0.6,0.8)-1e-10
ys= numpy.linspace(miny,maxy,201)
plot(ys,zvc(ys,E,pot=hp),color='k',zorder=10,lw=0.5)
plot(ys,-zvc(ys,E,pot=hp),color='k',zorder=11,lw=0.5);
Figure 13.17: Surface of section of the single chaotic orbit from Figure 13.16.
The chaotic orbit leaves many parts of the surface of section clear and we can identify the stable orbit families and the islands surrounding them directly. Of course, from this diagram we do not know whether the regions untouched by the chaotic orbits are regular or not, but we can investigate that by integrating orbits in each of these regions and determining their surface of section. Thus, the surface of section provides an easy way to investigate the orbital structure at a given energy.
.
.
.