Chapter 8, Section 8.12.1, Figure 8.22

The definitions of the following functions are taken from the PhD thesis of Malhotra (1988).

First we define the function ψ.

ψ[δ_] := N[ArcTan[Sqrt[(δ+1)^3-1]]];

Next we define the function Φ3.

Φ3[δ_] := N[(1/2)*(δ+1)*(Cos[ψ[δ]/3]+Sqrt[3]*Sin[ψ[δ]/3])^2];

The functions Φmin and Φmax determine the minimum and maximum values of Φ.  These are given in Eqs.(8.211) and (8.212).

Φmin[δ_] := N[3*(δ+1) - Φ3[δ] - 2*(2*Φ3[δ])^(1/4)];

Φmax[δ_] := N[3*(δ+1) - Φ3[δ] + 2*(2*Φ3[δ])^(1/4)];

The x values of the equilibrium points are the solutions of the cubic equation defined by xval.  Note that this is just Eq.(8.135) with the transformation given in Eq.(8.206).

xval[δ_]:=NSolve[x^3-3(δ+1)x-2 == 0, x];

The function Φval takes the real solutions given by xval and calculates the corresponding values of Φ.

Φval[δ_] := Sqrt[2*Max[Select[x /. xval[δ],(Im[#]==0)&]]]

Now define the Hamiltonian, h1, of the system.

h1[x_,y_,δ_]:=
   -(3/2)*(δ+1)*(x^2+y^2)+(x^2+y^2)^2/4-2*x;

The function xpts returns the values of x for y = 0 that lie on a particular curve of constant Hamiltonian.

xpts[hval_, δ_] := 
Select[x /. NSolve[h1[x,0,δ] == hval, x], (Im[#]==0)&]

The function ymax returns the maximum value of y for particular values of x, the Hamiltonian and δ.

ymax[x_, h_, δ_] :=
Module[{a,b,c},
a = (1/4);
b = (1/2)*x^2 - (3/2)*(δ+1);
c = (1/4)x^4 - (3/2)*(δ+1)*x^2 - 2*x - h;
Sqrt[(-b + Sqrt[b^2-4*a*c])/2/a]]

The function area calculates the area enclosed by a curve defined by the Hamiltonian value h for a particular δ.

area[h_,δ_] := 
Module[{xmin,xmax},
{xmin,xmax} = xpts[h,δ];
2*NIntegrate[y,{x,xmin,xmax},{y,0,ymax[x,h,δ]}]]

The function ΦValue returns the value of Φ for a particular δ and Hamiltonian h.

ΦValue[h_,δ_] := N[area[h,δ]/(2*Pi)]

The function h1Φ returns the Hamiltonian value taking input in terms of Φ, ϕ, and δ.

h1Φ[Φ_,ϕ_,δ_]:=N[-3(δ+1)*Φ + Φ^2 - 2*Sqrt[2*Φ]*Cos[ϕ]];

In the series of plots we are assuming that Φ, the action or alternatively the area enclosed by the trajectory, has a value of unity.  We start with δ = –6 and find a value of the Hamiltonian that gives a unit value for Φ.

test = ΦValue[22.9556,-6]

Now calculate the value of xmin and xmax for this value of the Hamiltonian.

{xmin, xmax} = xpts[22.9556,-6]

Now calculate the value of Φleft and Φright for these values.

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

Verify that this is indeed the correct value of the Hamiltonian.

h = N[h1Φ[Φright,0,-6]]

Plot the curve associated with this value of the Hamiltonian.  See Fig.8.22a.

p1 =
ContourPlot[h1[x,y,-6],{x,-5,5},{y,-5,5},
   ContourShading->False, AxesOrigin->{0,0},
   Contours->{h},PlotPoints->100]

Now go through the same procedure for δ = –1.

test = Chop[ΦValue[1.69027,-1]]

{xmin, xmax} = xpts[1.69027,-1]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,-1]]

Plot the curve associated with this value of the Hamiltonian.  See Fig.8.22b.

p2 =
ContourPlot[h1[x,y,-1],{x,-5,5},{y,-5,5},
   ContourShading->False,
   Contours->{h},PlotPoints->100]

Now go through the same procedure for δ = 0.

test = Chop[ΦValue[-2.38395,0]]

{xmin, xmax} = xpts[-2.38395,0]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,0]]

Here we also calculate the value of the Hamiltonian associated with the critical curve.

hcrit = h1Φ[Φmax[0],0,0]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.22c.

p3 =
ContourPlot[h1[x,y,0],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->100]

Now go through the same procedure for δ = 1.

test = Chop[ΦValue[-9.5931,1]]

{xmin, xmax} = xpts[-9.5931,1]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,1]]

hcrit = h1Φ[Φmax[1],0,1]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.22d.

p4 =
ContourPlot[h1[x,y,1],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->100]

Now go through the same procedure for δ = 3.

test = Chop[ΦValue[-37.12823,3]]

{xmin, xmax} = xpts[-37.12823,3]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,3]]

hcrit = h1Φ[Φmax[3],0,3]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.22e.

p5 =
ContourPlot[h1[x,y,3],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->100]