This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 2

Consider the following Hénon–Heiles-type potential (compared to the classical Hénon–Heiles-type potential, the additional (r^4) term guarantees confinement of the quantum mechanical problem)

V(x,y)=α ((((x^2)+(y^2)))^2)+β ((x^2)+(y^2))+γ ((x^3)-3x (y^2))=α (r^4)+β (r^2)+γ cos(3φ)(r^3).

For some values of the parameters α, β, and γ, and various initial conditions, calculate some of the classical trajectories. Visualize these trajectories and construct the corresponding Poincaré plots.

Also, calculate and visualize [✶], [✶] some of the energy eigenfunctions of the corresponding time-independent Schrödinger equation [✶], [✶]

(-((ℏ^2)/(2M))(Subscript[Δ, x,y])+V(x, y))(Subscript[Ψ, μ])(x, y)=(Subscript[ε, μ])(Subscript[Ψ, μ])(x, y)

by expanding the (Subscript[Ψ, μ])(x, y)=(∑_(m,n))(Subscript[c, m,n]^((μ)))(Subscript[ψ, n,m])(x, y) into the eigenfunctions (Subscript[ψ, n,m]) of the 2D harmonic oscillator

(-((ℏ^2)/(2M))(Subscript[Δ, r])+(M/2)ω ((|r|)^2))(Subscript[ψ, n,m])(r,φ)=(-((ℏ^2)/(2M))((1/r)(∂/(∂r))(r (∂/(∂r)))+(1/(r^2))((∂^2)/(∂(φ^2))))+(M/2)ω (r^2))(Subscript[ψ, n,m])(r,φ)
=(Subscript[ε, n,m])(Subscript[ψ, n,m])(r,φ)

for some values of the parameters α, β, and γ. The normalized eigenfunctions of the 2D harmonic oscillator in polar coordinates are (here (Subscript[|n,m⟩, r]) denotes the radial eigenfunctions with quantum numbers n and m and (Subscript[|m⟩, φ]) denotes the azimuthal eigenfunction with quantum number m).

⟨r,φ|n,m⟩=(Subscript[⟨φ|m⟩, φ])(Subscript[⟨r|n,m⟩, r])=(Subscript[ψ, n,m])(r,φ)=(Subscript[Φ, m])(φ)(Subscript[R, n,m])(r)
(Subscript[Φ, m])(φ)=(Subscript[⟨φ|m⟩, φ])=(1/(Sqrt[2π]))(e^(i m φ))
(Subscript[R, n,m])(r)=(Subscript[⟨r|n,m⟩, r])=((Sqrt[2])/L)(Sqrt[((n!)/((n+|m|)!))]) ((((r/L)))^(|m|)) (e^(-(1/2) (((r/L))^2))) (Subscript[L, n]^(|m|))(((((r/L)))^2))

⟨(Subscript[n, 2]),(Subscript[m, 2])|(Subscript[n, 1]),(Subscript[m, 1])⟩=(∫_0^(2π))(∫_0^∞)(Subscript[ψ, (Subscript[n, 1]),(Subscript[m, 1])])(r,φ)(OverBar[(Subscript[ψ, (Subscript[n, 2]),(Subscript[m, 2])])(r,φ)])r dr=1

and the eigenvalues are (Subscript[ε, n,m])=(2 n+|m|+1) ℏ ω [✶]. Here n=0,1,… and m=…,-2,-1,0,1,2,…. and L=(((ℏ/(M ω)))^(1/2)).

(* the generalized Hénon-Heiles potential *)
V[{x_, y_}, {α_, β_, γ_}] := α Sqrt[x^2 + y^2]^4 +
            β Sqrt[x^2 + y^2]^2 + γ (x^3 - 3 x y^2)

Show[GraphicsArray[
Block[{L1 = 6, L2 = 8, $DisplayFunction = Identity}, 
{(* contour plot of a concrete potential realization *)
 ContourPlot[Evaluate[V[{x, y}, {1, -6, -4}]], 
            {x, -L1, L1}, {y, -L1, L1}, PlotPoints -> 240, 
            Contours -> 60, ContourLines -> True, 
            (* inner contours more pronounced *)
            ContourStyle -> Table[{GrayLevel[x]}, {x, 0, 1, 1/60}],
            ColorFunction -> (GrayLevel[ArcCos[#]/(Pi/2)]&)],
 (* V(x,y) == 0 curves of the potential *)
 Show[Table[
 ContourPlot[Evaluate[V[{x, y}, {1, -6, γ}]], 
             {x, -L2, L2}, {y, -L2, L2},
             PlotPoints -> 240, Contours -> {0},
             ContourStyle -> {{Hue[0.78 (γ + 8)/16]}}, 
             ContourShading -> False], {γ, -8, 8}]]}]]]

Next, we visualize some of the classical trajectories. The function path calculates the path (trajectory) {x(t),y(t)} of a point particle in the potential V(x,y) with initial conditions {x(t),y(t)}={x0,y0}, {x'(t),y'(t)}={vx0,vy0}.

path[{x0_, y0_}, {vx0_, vy0_}, {α_, β_, γ_}, T_] :=
   NDSolve[(* equations of motion *)
   {x'[t] == vx[t], vx'[t] == -D[V[{x[t], y[t]}, {α, β, γ}], x[t]],
    y'[t] == vy[t], vy'[t] == -D[V[{x[t], y[t]}, {α, β, γ}], y[t]],
    x[0] == x0, y[0] == y0, vx[0] == vx0, vy[0] == vy0},
   {x, y, vx, vy}, {t, 0, T}, MaxSteps -> 10^6]

We visualize the path using the function pathPlots. For a given path pathSol (returned from the function path), this function generates the four parametric plots {x'(t),y'(t)}, {x(t),x'(t)}, {y(t),y'(t)}, and {x'(t),y'(t)},

pathPlots[pathSol_] :=
Module[{X, Y, VX, VY},
  (* extract value lists for positions and velocities *)
  {X, Y, VX, VY} = 
      Transpose[Partition[#[[2, 4, 2]], 2]][[1]]& /@ pathSol[[1]];
Show[GraphicsArray[
  ListPlot[Transpose[#], PlotJoined -> True,
         PlotStyle -> {Thickness[0.002], GrayLevel[0]},
         PlotRange -> All, (* for uniformity *) AspectRatio -> 1,
         DisplayFunction -> Identity, FrameTicks -> False,
         Frame -> True, PlotLabel -> #2, Axes -> False]& @@@ 
         (* the four pairs of functions for the parametric plots *)
         {{{X, Y}, {x, y}}, {{X, VX}, {x, Subscript[v, x]}},
          {{Y, VY}, {y, Subscript[v, y]}},
          {{VX, VY}, {Subscript[v, y], Subscript[v, xy]}}}]]]

Here are some of the trajectories and their phase plots shown.

vβδData =
(* parameter sets for various types of classical motion *)
{{{-1475/496,    66/041}, {   19/003, - 326/221}},
 {{ 5043/515,   906/235}, {-1200/157, -3576/569}},
 {{ 2923/321,  1517/404}, { 1690/317, -  23/021}},
 {{- 524/147,  2272/567}, {-  58/043,     9/035}},
 {{ 6097/823, -2230/233}, {- 648/271, -   1/107}},
 {{ 2305/634, - 221/133}, { 5874/959,  1287/458}},
 {{  753/106,   349/150}, {- 665/624, -5509/732}}};

pathPlots[path[{0, 0}, #2, {1, Sequence @@ #1}, 50]]& @@@ vβδData

Here are 3D plots of some of the trajectories in x,y,(Subscript[v, x])-space.

pathPlot3D[pathSol_] :=
Module[{X, Y, VX, VY, trajectory, λ},
  (* extract value lists for positions and velocities *)
  {X, Y, VX, VY} = 
      Transpose[Partition[#[[2, 4, 2]], 2]][[1]]& /@ pathSol[[1]];
  (* path in x,y,vx-space *) 
  xyvxPath = Transpose[{X, Y, VX}]; λ = Length[xyvxPath];    
Show[Graphics3D[{Thickness[0.002],
       (* color path from red to blue *)
       MapIndexed[{Hue[0.78 #2[[1]]/λ], Line[#1]}&, 
                   Partition[xyvxPath, 2, 1]]}],
     BoxRatios -> {1, 1, 1}, Axes -> True,
     AxesLabel -> {"x", "y", Subscript[v, x]}]]

Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
pathPlot3D[path[{0, 0}, #2, {1, Sequence @@ #1}, 100]]& @@@ 
                                    vβδData[[{1, 6, 7}]]]]]

Some of the trajectories look regular, and some look chaotic [✶], [✶], [✶], [✶]. Displaying their Poincaré plots confirms this [✶], [✶]. For a given pathSol returned from the function path, the function poincareSection generates the Poincaré plot of the variables {ξ,η} for ζ(t)=0.

poincareSection[pathSol_, {ξ_, η_}, ζ_] :=
Module[{ζL, tL, tζL, tζLPairs, ts, points}, 
 (* find times where ζ vanishes *)
 ζL = Transpose[Partition[(ζ /. pathSol[[1]])[[4, 2]], 2]][[1]];
 tL = (ζ /. pathSol[[1]])[[3, 1]];
 tζL = Transpose[Partition[#, 2, 1]& /@ {tL, ζL}];
 tζLPairs = Select[tζL, (#[[2, 1]] #[[2, 2]] < 0)&];
 (* use linear interpolation *)
 ts = (#1[[1]] + #2[[1]]/(#2[[1]] - #2[[2]])*
       (#1[[2]] - #1[[1]]))& @@@ tζLPairs;
 (* {ξ, η} at the times ts *)
 points = ({ξ[#], η[#]}& /@ ts) /. pathSol[[1]];
 (* form Poincare plot graphics *)
 Graphics[{PointSize[0.003], Point /@ points},
       PlotRange -> All, AspectRatio -> 1, Frame -> True, 
       FrameTicks -> False, PlotLabel -> {ξ, η}, Axes -> False]]

For a given path, the function poincarePlots shows the four Poincaré plots {x,y} at x'(t)=0, {x,y} at y'(t)=0, {x',y'} at x(t)=0, and {x',y'} at x(t)=0.

poincarePlots[pathSol_] := 
Show[GraphicsArray[poincareSection[pathSol, #1, #2]& @@@
{{{x, y}, vx}, {{x, y}, vy}, {{vx, vy}, x}, {{vx, vy}, y}}]]

Here are the Poincaré plots corresponding to the above parametric and phase plots. To obtain enough intersections of the coordinate planes with the trajectory, we now propagate to the much larger time t=5000.

poincarePlots[path[{0, 0}, #2, 
                   {1, Sequence @@ #1}, 5000]]& @@@ vβδData

(For closed-form solutions of the classical Hénon–Heiles potential, see [✶], [✶], [✶].)

We continue with the quantum mechanical treatment of the generalized Hénon–Heiles potential. To calculate the quantum mechanical eigenfunctions, we have to diagonalize the Hamiltonian corresponding to V(x,y) in a truncated 2D harmonic oscillator basis. We write the matrix elements as a sum of the harmonic oscillator diagonal terms and perturbation terms. (This avoids the need for the explicit calculation of the matrix elements of the kinetic part in exchange for the anyway-to-be-calculated matrix element of the (r^2) term.) So we have (using the notations (Subscript[H, 0]) and (Subscript[V, 0]) for the Hamiltonian and the potential of the 2D harmonic oscillator) the following expression for a generic matrix element:

⟨(Subscript[n, 1]),(Subscript[m, 1])|-((ℏ^2)/(2M))Δ+V|(Subscript[n, 2]),(Subscript[m, 2])⟩=⟨(Subscript[n, 1]),(Subscript[m, 1])|(Subscript[H, 0])+V-(Subscript[V, 0])|(Subscript[n, 2]),(Subscript[m, 2])⟩
=(2(Subscript[n, 1])+(Subscript[m, 1])+1)(Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])])(Subscript[δ, (Subscript[n, 1]),(Subscript[n, 2])])ℏ ω+⟨(Subscript[n, 1]),(Subscript[m, 1])|V-((M )/2)(ω^2)(r^2)|(Subscript[n, 2]),(Subscript[m, 2])⟩=
(2(Subscript[n, 1])+m+1)(Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])])(Subscript[δ, (Subscript[n, 1]),(Subscript[n, 2])])ℏ ω+ ⟨(Subscript[n, 1]),(Subscript[m, 1])|α (r^4)+β (r^2)-((M )/2)(ω^2)(r^2)|(Subscript[n, 2]),(Subscript[m, 2])⟩+
⟨(Subscript[n, 1]),(Subscript[m, 1])|γ cos(3φ)(r^3)|(Subscript[n, 2]),(Subscript[m, 2])⟩.

The φ-independent potential term is diagonal in the angular quantum numbers (Subscript[m, 1]) and (Subscript[m, 2]).

⟨(Subscript[n, 1]),(Subscript[m, 1])|α (r^4)+β (r^2)-((M )/2)(ω^2)(r^2)|(Subscript[n, 2]),(Subscript[m, 2])⟩= (Subscript[(Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])]) , r])(Subscript[⟨(Subscript[n, 1]),(Subscript[m, 1])|α (r^4)+β (r^2)-((M )/2)(ω^2)(r^2)|(Subscript[n, 2]),(Subscript[m, 2])⟩, r]).

The φ-dependent potential term factors into an angular and an radial part.

⟨(Subscript[n, 1]),(Subscript[m, 1])|cos(3φ)(r^3)|(Subscript[n, 2]),(Subscript[m, 2])⟩ =(Subscript["", φ])(Subscript[(Subscript[⟨(Subscript[m, 1])|cos(3φ)|(Subscript[m, 2])⟩, φ])   , r])(Subscript[⟨(Subscript[n, 1]),(Subscript[m, 1])|γ (r^3)|(Subscript[n, 2]),(Subscript[m, 2])⟩, r]).

We anticipate (from the corresponding 1D examples of Chapter 2) needing a few hundred angular and a few hundred radial basis states of the harmonic oscillator. This means we will have in the order of a few (10^4) product basis states and the resulting Hill matrix will have in the order (10^9) matrix elements. Calculating (or storing and manipulating) (10^9) matrix elements individually is, of course, hopeless. However, being an exercise, we anticipate that for a given set of quantum numbers (Subscript[n, 1]), (Subscript[m, 1]), only finitely many states |(Subscript[n, 2]),(Subscript[m, 2])⟩ will have nontrivial matrix elements with V(x,y). These nonvanishing matrix elements we will use to create and diagonalize a sparse matrix. Calculating explicit forms of these matrix elements is the next step.

Instead of using the complex-valued eigenfunctions exp(±i m φ), we will use the real-valued ones cos (m φ) and sin(m φ). While this slightly complicates some of the following calculations, we gain two advantages: a) Because the potential V(x, y) is mirror symmetric with respect to the x-axis, so will be the eigenfunctions. This means in the eigenfunction expansions we can use a smaller matrix size to take a fixed number of basis states into account. b) All arithmetic needed will only operate on real numbers, which is considerably faster than complex number calculations. An additional small benefit of the division in even and odd states is the possibility to optimize the graphics constructions later when we visualize the eigenfunctions.

The basic angular orthogonality integrals are

(∫_0^(2π))cos((Subscript[m, 1])φ)cos((Subscript[m, 2])φ)dφ=(1+(Subscript[δ, (Subscript[m, 1]),0]))(Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])])π     for (Subscript[m, 1]), (Subscript[m, 2])=0,1,2,…

(∫_0^(2π))sin((Subscript[m, 1])φ)sin((Subscript[m, 2])φ)dφ=(Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])])π     for (Subscript[m, 1]), (Subscript[m, 2])=1,2,…

and the basic radial orthogonality integral is

(∫_0^∞)(t^λ) (e^(-t)) (Subscript[L, (Subscript[n, 1])]^λ)(t) (Subscript[L, (Subscript[n, 2])]^λ)(t)dt=((((Subscript[n, 1])+λ+1)! )/((Subscript[n, 1])!))(Subscript[δ, (Subscript[n, 1]),(Subscript[n, 2])]).

To reduce products of powers of r with the Laguerre polynomials to this standard form, we make use of the following two recursion relations for associated Laguerre polynomials (Subscript[L, n]^λ)(z) [✶]:

z (Subscript[L, n]^λ)(z)=-(n+λ) (Subscript[L, n-1]^λ)(z)+(2 n+λ+1) (Subscript[L, n]^λ)(z)-(n+1) (Subscript[L, n+1]^λ)(z)
(Subscript[L, n]^(λ+1))(z)=(λ+n) (Subscript[L, n-1]^λ)(z)+(z-n) (Subscript[L, n]^λ)(z).

Off[General::spell1]; Off[General::spell1];

Here are the even and odd eigenfunctions of the 2D harmonic oscillator. Because later we will set ℏ=M=1, we do not include these variables as arguments in the left-hand sides of the corresponding function definitions for R, ψEven, and ψOdd. We will set a value for ω later globally (to avoid a cluttered notation we also do not use ω as an explicit argument).

(* natural length scale *)
L = Sqrt[ℏ/(M ω)];

(* angular momentum eigenfunctions *)
ΦOdd[m_, φ_] = 1/Sqrt[Pi] Sin[m φ];
(* normalization for even states with m = 0 is special *)
ΦEven[m_, φ_] = (1 - (1 - 1/Sqrt[2])  KroneckerDelta[m]) *
                1/Sqrt[Pi] Cos[m φ];

(* radial momentum eigenfunctions *)
R[n_, m_, r_] := Sqrt[2]/L Sqrt[n!/(n + m)!] (r/L)^m*
                 Exp[-(r/L)^2/2] LaguerreL[n, m, (r/L)^2];
    
(* 2D eigenfunctions *)
ψOdd[n_, m_, r_, φ_] := ΦOdd[m, φ] R[n, m, r];
ψEven[n_, m_, r_, φ_] := ΦEven[m, φ] R[n, m, r];

Without loss of generality, we can later set ℏ=M=1. While the eigenfunctions are complete for every values of ω, we will not set ω to a specific value, but rather keep it as a parameter that will allow to radially expand and contract the basis.

Here is a quick check of the correctness and normalization of the eigenfunctions.

SchrödingerOperatorHO2D[ψ_, {r_, φ_}] :=
-ℏ^2/(2 M) (1/r D[r D[ψ, r], r] + 1/r^2 D[ψ, φ, φ]) + 1/2 M ω^2 r^2 ψ

Block[{psi},
Table[psi = ψEven[n, m, r, φ]; 
      (* check if Schrödinger equation is fulfilled *)
      Expand[SchrödingerOperatorHO2D[psi, {r, φ}] - 
      (2n + Abs[m] + 1) ℏ ω psi] == 0, {n, 0, 2}, {m, 0, 2}]] //
                                                Flatten // Union

Block[{ℏ = 1, M = 1, ω = 1},
Table[(* check normalization of eigenfunctions *)
      Integrate[ψEven[n1, m1, r, φ] ψEven[n2, m2, r, φ] r, 
                {φ, 0, 2Pi}, {r, 0, Infinity}] ==
      KroneckerDelta[n1, n2] KroneckerDelta[m1, m2],
      {n1, 0, 2}, {m1, 0, 2}, {n2, 0, n1}, {m2, 0, m1}]] //
                                               Flatten // Union

As a side step, let us have a quick look at the radial eigenfunctions. For increasing n and m, the eigenfunctions extend further radially outwards.

Block[{M = 1, ℏ = 1, ω = 1},
 Show[Graphics3D[
      Table[{Hue[0.8 m/(3 n)], 
             Line[Table[Evaluate[{r, n, R[n, m, r]^2}], 
                        {r, 0, 10, 1/10}]]},
            {n, 0, 8}, {m, (* odd states *) 1, 3 n}]],
            BoxRatios -> {2.4, 1, 1}, PlotRange -> All,
            Axes -> True, AxesLabel -> {"r", "n", ψ^2}]]

The next graphic shows the position of the maxima of R[n, m, r] for a fixed n=12 as a function of m.

radialMaxima[n_, m_, rMax_] :=
Block[{M = 1, ℏ = 1, ω = 1},
Point[{#[[2, 1]], Log[#[[2, 2]]]}]& /@ 
Select[Partition[Plot[Evaluate[R[n, m, r]^2], {r, 0, rMax},
                      PlotPoints -> 20,
                      DisplayFunction -> Identity][[1, 1, 1, 1]], 
                 3, 1], (#[[1, 2]] < #[[2, 2]] > #[[3, 2]])&]]

Block[{n = 12},
Show[Graphics[
      Table[{Hue[0.8 m/(8n)], radialMaxima[n, m, 16]},
            {m, 1, 8n, 4}]], PlotRange -> {0, -5}, Axes -> True]]

We have to decide which harmonic oscillator eigenfunctions to include in the truncated basis. We have a doubly infinite set of eigenfunctions. The above graphics showed that for an increasing n, the eigenfunctions expand increasingly radially outwards. The generalized Hénon–Heiles potential is bounded in asymptotically increasing in every direction due to the (r^4) term. This means, all the energy eigenfunctions will be well localized in the corresponding classically allowed region. So we fix an energy εo, and take all eigenfunctions having an energy less than or equal to εo into account. We sort these selected eigenfunctions by energy, and in case of degeneracy, by increasing m. The functions numberEvenBasisFunctions and numberOddBasisFunctions carry out this selection and ordering process. They create the globally visible functions efNumber and nmNumbers and the variable dim. For a given pair of quantum numbers n, m, the function efNumber returns the position of this state in the basis. For a given position 1, 2, …, in the basis state, the function nmNumbers returns the pair of quantum numbers n, m. dim is the dimension of the basis. (For a cleaner, “ production code” quality programming one could export pure functions with a dispatched rule table in the body for efNumber, to keep all parts modular and well separated. But for the purpose of visualizing some eigenstates, which is the sole goal of this exercise, the current approach is appropriate.)

numberEvenBasisFunctions[εo_] :=
Module[{energyOrderedQuantumNumbers, nmPairs},
 (* order pairs of quantum numbers {n, m} by energy *)
 energyOrderedQuantumNumbers = 
  {#[[1, 2]], Sort[First /@ #]}& /@ Split[Sort[Select[Flatten[
  Table[{{n, m}, 2 n + m + 1}, 
        {n, 0, εo/2}, {m, 0, εo}], 1], #[[-1]] <= εo&],
        #1[[-1]] < #2[[-1]]&], #1[[-1]] == #2[[-1]]&];
 (* export the functions efNumber and nmNumbers the number dim *)
 Clear[efNumber, nmNumbers]; 
 (efNumber[#1] = #2; nmNumbers[#2] = #1)& @@@
 MapIndexed[{#1, #2[[1]]}&, nmPairs =
            Flatten[Last /@ energyOrderedQuantumNumbers, 1]];
 dim = Length[nmPairs]]

(* no m == 0 states in the odd functions *)
DownValues[numberOddBasisFunctions] =
 DownValues[numberEvenBasisFunctions] /.
  {numberEvenBasisFunctions -> numberOddBasisFunctions,
   HoldPattern[{m, 0, εo}] :> {m, 1, εo}};

Here is an example of the ordered n,m quantum numbers up to energy 12ℏω. The odd basis misses the states with m=0.

numberEvenBasisFunctions[12];
nmNumbers /@ Range[dim]

numberOddBasisFunctions[12];
nmNumbers /@ Range[dim]

The number of basis states grows quadratically with the energy εo. With energies in the few hundreds, we have a few ten thousand basis states.

(* count even and odd  basis sizes as a function of εo *)
Function[εo, {εo, numberEvenBasisFunctions[εo]; dim,
                  numberOddBasisFunctions[εo]; dim}] /@
                      {1, 5, 10, 50, 100, 200, 250, 500}

Because of the virial theorem, the average squared distance of the state ⟨r,φ|n,m⟩ is proportional to the energy (Subscript[ε, n,m]). This means that the ordering of the states carried out by the functions numberEvenBasisFunctions and numberEvenBasisFunctions can model states efficiently with an approximately known spatial extension.

rSquareMean[{n_, m_}] := 
Block[{M = 1, ℏ = 1, ω = 1},
      (* mean square spatial extension *)
      Integrate[r^2 R[n, m, r]^2 r, {r, 0, Infinity}]]

numberEvenBasisFunctions[12];
rSquareMean /@ (nmNumbers /@ Range[dim])

We use the sorted-by-energy states for making one more graphic—the integrated local density of states as a function of the radius and the energy up to (Subscript[ε, max])=36. The left graphic shows the integrated (over the energy) local density of states. The emergence of the classical limit is clearly visible. The right graphic shows the oscillations (due to the oscillating nature of the Laguerre polynomials) of the density of states better by subtracting the nonoscillating Thomas–Fermi contribution (∼(ε-V(r))θ(ε-V(r))) [✶].

Block[{εMax = 36, rMax = 10, (* avoid 0^0 *) rMin = 10^-6, 
       ppr = 160, ℏ = 1, M = 1, ω = 1,
       nmGroups, ψData, ldos, polysH, polysV, opts},
 (* the eigenfunctions up to energy εMax *)
 numberEvenBasisFunctions[εMax];
 nmGroups = Map[Last, Split[{2#1 + #2 + 1, {##}}& @@@ 
                (nmNumbers /@ Range[dim]), #1[[1]] === #2[[1]]&], {2}];
 (* summed squared eigenfunctions times r *)
 ψData = FoldList[Plus, Table[0, {r, rMin, rMax, (rMax - rMin)/ppr}],  
          (* sum eigenfunctions for a fixed energy *) Plus @@@ 
           Apply[Table[If[#2 === 0, 1/2, 1] R[#1, #2, N[r]]^2, 
                 {r, rMin, rMax, (rMax - rMin)/ppr}]&, nmGroups, {2}]];
 (* local density of states *)
 ldos = MapIndexed[{rMin + (#2[[2]] - 1)/ppr (rMax - rMin), 
                    #2[[1]] - 1, #1}&, ψData, {2}]; LD = ldos;
 polysH = (* horizontal polygons *)
 Table[Polygon[{ldos[[k, i]], ldos[[k, i + 1]],
                ldos[[k, i + 1]] + {0, 1, 0},
                ldos[[k, i]] + {0, 1, 0}}], 
       {k, Length[ldos]}, {i, Length[ldos[[k]]] - 1}];
 polysV = (* vertical polygons *)
 Table[Polygon[{ldos[[k, i]], ldos[[k, i + 1]],
                ldos[[k + 1, i + 1]], ldos[[k + 1, i]]}], 
       {k, Length[ldos] - 1}, {i, Length[ldos[[k]]] - 1}];
 (* graphics of the local density of states *)
 opts = Sequence[PlotRange -> All, Axes -> True, 
                 BoxRatios -> {2, 2, 1}];
 Show[GraphicsArray[
 {(* total integrated density of states *)
  Graphics3D[{EdgeForm[], polysH, polysV}, 
             AxesLabel -> {"r", "ε", "𝒟"}, opts],
  (* oscillating part of the total integrated density of states *)
  Graphics3D[{EdgeForm[], polysH, polysV} /. 
   Polygon[l_] :> Polygon[Function[{r, ε, 𝒟}, 
    (* with subtracted Thomas-Fermi contribution *)
    {r, ε, 𝒟 - 1/2 (ε - r^2/2) UnitStep[ε - r^2/2]}] @@@ l], 
             AxesLabel -> {"r", "ε", "δ𝒟"}, opts]}]]]

Next, we calculate the actual matrix elements. We start by calculating the matrix elements arising from the azimuthal eigenfunctions. For the φ-independent parts of the potential the matrix elements are obviously (Subscript[δ, (Subscript[m, 1]),(Subscript[m, 2])]). For the φ-dependent part, the function φMatrixElement calculates matrix elements in the form of Kronecker delta functions.

φMatrixElement[{Φ1_, Φ2_}, {m1_, m2_}, φPart_, φ_] := 
   TrigReduce[Φ1 φPart Φ2] /. (* emulate integration *)
   (* argument of Cos must vanish for nonvanishing integral *)
   Cos[x_] :> With[{𝓂2 = m2 /. Solve[x == 0, m2][[1]]},
                    2 Pi KroneckerDelta[m2, 𝓂2]] // Expand

For the odd-odd matrix elements, we see nonvanishing terms of the form (Subscript[δ, (Subscript[m, 1]),(Subscript[m, 1])±3]) and (Subscript[δ, 1,2]) and, because of symmetry (Subscript[δ, 2, 1]).

φMatrixElement[{ΦOdd[m1, φ], ΦOdd[m2, φ]}, {m1, m2}, Cos[3 φ], φ]

For a specific value of (Subscript[m, 1]), the function MEOddφList returns a list of the values of (Subscript[m, 2]) so that we have a nonvanishing matrix element and its value.

MEOddφList[m1_?Positive] :=
Cases[{{3 - m1, -1/2}, {m1 + 3, 1/2}, {m1 - 3, 1/2}}, {_?Positive, _}]

The even-even matrix elements have a similar structure. Because (Subscript[Φ, 0]^((e)))(φ) is now nonvanishing, we have also terms of the form (Subscript[δ, 0,3]) and (Subscript[δ, 3,0]) (which, because of the different normalization for m=0, are slightly different from (Subscript[δ, (Subscript[m, 1]),(Subscript[m, 1])±3]).

φMatrixElement[{ΦEven[m1, φ], ΦEven[m2, φ]}, {m1, m2}, 
               Cos[3 φ], φ] // Simplify

For a specific value of (Subscript[m, 1]), the function MEEvenφList returns a list of the values of (Subscript[m, 2]) so that we have a nonvanishing matrix element and its value.

MEEvenφList[m1_?NonNegative] :=
Cases[{If[m1 == 0, {3, 1/Sqrt[2]}, {}],
       If[m1 == 3, Sequence @@ {{0, 1/Sqrt[2]}, {m1 + 3, 1/2}}, {}],
       If[m1 != 0 && m1 != 3, 
          Sequence @@ {{3 - m1, 1/2}, {m1 + 3, 1/2}, {m1 - 3, 1/2}}, {}]},
      {_?NonNegative, _}]

The matrix elements between the even and the odd states vanish for the cos(3φ) potential part.

TrigReduce[ΦEven[m1, φ] Cos[3 φ] ΦOdd[m2, φ]] /. _Sin -> 0

Next, we calculate the integrals arising from the radial part of the matrix elements of the potential. For the φ-independent potential part (meaning (Subscript[m, 1])=(Subscript[m, 2])=m), we have the two integrals (∫_0^∞)(r^ε)(Subscript[R, (Subscript[n, 1]),m])(r)  (Subscript[R, (Subscript[n, 2]),m])(r) r dr with ε=2 and ε=4. After a change of variables r->L (z^(1/2)) these reduce to the following type integral [✶], [✶]

(∫_0^(∞ ))(r^ε)(Subscript[R, (Subscript[n, 1]),m])(r)  (Subscript[R, (Subscript[n, 2]),m])(r) r dr=(L^ε)(Sqrt[(((Subscript[n, 1])!(Subscript[n, 2])!)/(((Subscript[n, 1])+|m|)!((Subscript[n, 2])+|m|)!))]) (∫_0^(∞ ))(z^(ε/2))(e^(- (z^2))) (z^m)(Subscript[L, (Subscript[n, 1])]^m)(z)(Subscript[L, (Subscript[n, 2])]^m)(z) dz.

For the φ-dependent potential part (meaning (Subscript[m, 2])=(Subscript[m, 1])±3 or (Subscript[m, 2])=(Subscript[m, 1])±1), we have the two integral (∫_0^∞)(r^3)(Subscript[R, (Subscript[n, 1]),m])(r)  (Subscript[R, (Subscript[n, 2]),m±μ])(r) r dr with μ=3 and μ=1 (for some special states (Subscript[n, 1,2])=1,2). Carrying out again the change of variables r->L (z^(1/2)) these integrals reduces to the following (shown are only the cases of positive μ, the cases of negative μ follow from a simple shift in m).

(∫_0^(∞ ))(r^3)(Subscript[R, (Subscript[n, 1]),m])(r)  (Subscript[R, (Subscript[n, 2]),m+μ])(r) r dr=(L^((3+μ)/2))(Sqrt[(((Subscript[n, 1])!(Subscript[n, 2])!)/(((Subscript[n, 1])+m)!((Subscript[n, 2])+m+μ)!))]) (∫_0^(∞ ))(z^((3+μ)/2))(e^(- (z^2))) (z^m)(Subscript[L, (Subscript[n, 1])]^m)(z)(Subscript[L, (Subscript[n, 2])]^(m+μ))(z) dz

Using the above recursion relations for z (Subscript[L, n]^λ)(z) and (Subscript[L, n]^(λ+1))(z) repeatedly, we now reduce the last two integrals to standard integrals of the orthogonality type (∫_0^∞)(t^λ) (e^(-t)) (Subscript[L, n]^λ)(t) (Subscript[L, p]^λ)(t)dt∼(Subscript[δ, n,p]).

The two functions LaguerreReduceRule1 and LaguerreReduceRule2 apply the two recursion relations repeatedly.

(* remove power prefactor by shifting Laguerre order *)
LaguerreReduceRule1 = z^(μ_ + k_Integer) LaguerreL[n_, μ_, z] :> 
  -(n + μ) z^(μ + k - 1) LaguerreL[n - 1, μ, z] +  
   (1 + 2 n + μ) z^(μ + k - 1) LaguerreL[n, μ, z] -
   (n + 1) z^(μ + k - 1) LaguerreL[n + 1, μ, z];

(* remove shifted parameter by shifting Laguerre order *)
LaguerreReduceRule2 = LaguerreL[n_, μ_ + i_Integer, z] :> 
   (n + (μ + i - 1))/z LaguerreL[n - 1, μ + i - 1, z] + 
   (z - n)/z LaguerreL[n, μ + i - 1, z];

Here is an example showing the two functions at work. We canonicalize the expression (z^(m+4))(Subscript[L, n]^(m+2)). The summands after expansion all have the desired form (z^m)(Subscript[L, ν]^m)(z).

Collect[Expand[z^(m + 4) LaguerreL[n, m + 2, z] //. 
                      LaguerreReduceRule2] //. LaguerreReduceRule1,
        _LaguerreL, Factor]

The function LaguerreOverlapIntegral forms the actual overlap integral (∫_0^∞)zrPart (e^(-z)) (Subscript[L, (Subscript[n, 1])]^m)(z) (Subscript[L, (Subscript[n, 2])]^(m+i))(z)dz.

LaguerreOverlapIntegral[{{n1_, m_}, {n2_, m_ + i_.}}, zrPart_] := 
Collect[(* apply the two recursion relations *)
        Expand[Expand[zrPart LaguerreL[n2, m + i, z] //. 
                      LaguerreReduceRule2] //. LaguerreReduceRule1],
        _LaguerreL, Factor] /. 
        (* use basic orthogonality integral *)
        z^μ_ LaguerreL[n_, μ_, z] :> 
         (n1 + μ)!/n1! KroneckerDelta[n1, n] /. μ -> m

The function δRule finally expresses (Subscript[n, 2]) as a function of (Subscript[n, 1]) in the Kronecker symbols. This allows, for a given (Subscript[n, 1]), to explicitly list all nonvanishing matrix elements.

δRule[expr_] := (List @@ Collect[Expand[expr], 
                                 _KroneckerDelta, Factor]) /.
 a_ KroneckerDelta[n1_, ν_] :> ({n2, a} /. Solve[n1 == ν, n2][[1]]);

Now, we have all tools together to efficiently evaluate the radial matrix elements of the potential. The function MEr calculates the radial integrals belonging to the α (r^4)+β (r^2) potential part.

MEr[n1_, n2_, m_, rPart_] := 
Sqrt[n1!/(n1 + m)!] Sqrt[n2!/(n2 + m)!]*
 L^(Exponent[rPart, r]) LaguerreOverlapIntegral[{{n1, m}, {n2, m}}, 
                                         z^m z^(Exponent[rPart, r]/2)]

For later use in the construction of the sparse global Hill matrix, we write all nonvanishing matrix elements in the form {(Subscript[n, 2])((Subscript[n, 1])), value((Subscript[n, 1]))}. These are the corresponding matrix elements for the quadratic potential term.

MErList1[{n1_, m_}] = 
   δRule[(-M/2 ω^2 + β) MEr[n1, n2, m, r^2]] // PowerExpand

Here are the matrix elements for the quartic potential term.

MErList2[{n1_, m_}] = δRule[α MEr[n1, n2, m, r^4]] // PowerExpand

The function MEr3 calculates the nonvanishing radial overlap integrals for the potential part γ (r^3) and (Subscript[m, 2])=(Subscript[m, 1])±3.

MEr3[n1_, n2_, m_, rPart_] := 
Sqrt[n1!/(n1 + m)!] Sqrt[n2!/(n2 + m + 3)!]*
 L^(Exponent[rPart, r] - 3)*
    LaguerreOverlapIntegral[{{n1, m}, {n2, m + 3}}, 
                            z^m z^(Exponent[rPart, r]/2)]

Here are these nonvanishing matrix elements.

MErListP3[{n1_, m_}] = δRule[γ MEr3[n1, n2, m, r^6] ] // PowerExpand

After some changes of variables, the case (Subscript[m, 2])=(Subscript[m, 1])-3 reduces to the case (Subscript[m, 2])=(Subscript[m, 1])+3.

MErListM3[{n1_, m_}] = δRule[γ MEr3[n1, n2, m, r^6] /.
               {n1 -> n2, n2 -> n1, m -> m - 3}] // PowerExpand

The function MEr12 finally calculates the nonvanishing radial overlap integrals for the potential part γ (r^3) and (Subscript[m, 2])=(Subscript[m, 1])±1.

MEr12[n1_, n2_, m_, rPart_] := 
Sqrt[n1!/(n1 + m)!] Sqrt[n2!/(n2 + m + 1)!]*
 L^(Exponent[rPart, r] - 1)*
    LaguerreOverlapIntegral[{{n1, m}, {n2, m + 1}}, 
                                         z^m z^(Exponent[rPart, r]/2)]

Here are these nonvanishing matrix elements.

MErListP12[{n1_, m_}] = δRule[γ MEr12[n1, n2, m, r^4] ] // PowerExpand

MErListM12[{n1_, m_}] = δRule[γ MEr12[n1, n2, m, r^4] /.
               {n1 -> n2, n2 -> n1, m -> m - 1}] // PowerExpand

Now, we join the last results, and incorporate the values of the angular overlap integrals. For given quantum numbers (Subscript[n, 1]) and (Subscript[m, 1]), the functions nonvanishingOddMatrixElements and nonvanishingOddMatrixElements return a list of all nonvanishing matrix elements of the form {{(Subscript[n, 2])((Subscript[n, 1])), (Subscript[m, 2])((Subscript[n, 1]))}, matrixElement}. We now set ℏ=M=1, but keep ω still as a symbolic parameter. Because the even and odd case have similar structures, we define an intermediate function nonvanishingMatrixElements.

(* treat even and odd case uniformly *)
nonvanishingMatrixElements[{n1_, m1_}, {𝒻1_, 𝒻2_, ℊ_, 𝒽_}, 
                           {α_, β_, γ_}, ω_] =
 (* incorporate φ-matrix element values as prefactors *)
 Join[(* return in the form {n2(n1, m1), m2(m1)}, matrixElement} *)
 {(* harmonic oscillator part *)
  {{n1, m1}, (2 n1 + m1 + 1) ℏ ω}},
   (* r-dependent potential parts *)
   If[#1 >= 0, {{#1, m1}, #2}, {}]& @@@ MErList1[{n1, m1}], 
   If[#1 >= 0, {{#1, m1}, #2}, {}]& @@@ MErList2[{n1, m1}], 
   (* cos(3φ) potential part *)
   If[#1 >= 0, {{#1, m1 + 3}, 𝒻1 #2}, {}]& @@@ MErListP3[{n1, m1}], 
   If[#1 >= 0 && (m1 ~ 𝒽 ~ 3), {{#1, m1 - 3}, 𝒻2 #2}, {}]& @@@ 
                                              MErListM3[{n1, m1}],
   If[m1 == 1 && #1 >= 0, {{#1, 2}, ℊ #2}, {}]& @@@ 
                                              MErListP12[{n1, m1}],
   If[m1 == 2 && #1 >= 0, {{#1, 1}, ℊ #2}, {}]& @@@ 
                                              MErListM12[{n1, m1}]] /. 
                                                    {M -> 1, ℏ -> 1};

The differences between the even and the odd cases are the inclusion/exclusion of the m=0 case and different signs of the angular momentum matrix elements.

(* even case *)
nonvanishingEvenMatrixElements[{n1_, m1_}, {α_, β_, γ_}, ω_] =
 nonvanishingMatrixElements[{n1, m1}, (* special rule for m == 0 *)
             {If[m1 == 0, 1/Sqrt[2], 1/2], 
              If[m1 == 3, 1/Sqrt[2], 1/2], 1/2,
              GreaterEqual}, {α, β, γ}, ω];

(* odd case *)
nonvanishingOddMatrixElements[{n1_, m1_}, {α_, β_, γ_}, ω_] =
 nonvanishingMatrixElements[
             {n1, m1}, {1/2, 1/2, -1/2, Greater}, {α, β, γ}, ω];

Here is a shortened form of the function nonvanishingEvenMatrixElements.

nonvanishingEvenMatrixElements[{n1, m1}, {α, β, γ}, ω] //
                                            Short[#, 12]&

For explicitly given values of (Subscript[n, 1]) and (Subscript[m, 1]), the test in nonvanishingMatrixElements can completely be carried out. As a result, we can add all nonvanishing matrix element values belonging to identical {(Subscript[n, 2]),(Subscript[m, 2])} pairs. The function uniteElements does this. The two function allNonvanishingOddMatrixElements and allNonvanishingOddMatrixElements are the final functions used to construct the Hill matrices.

uniteElements[l_] := (* sum right-hand side of rules *)
If[Length[#] === 1, Rule @@ First[#], 
   #[[1, 1]] -> (Plus @@ (Last /@ #))]& /@ 
     Split[Sort[DeleteCases[l, {}]], #1[[1]] === #2[[1]]&]

allNonvanishingEvenMatrixElements[{n1_Integer, m1_Integer}, 
                                  {α_, β_, γ_}, ω_] := 
uniteElements @ nonvanishingEvenMatrixElements[{n1, m1}, {α, β, γ}, ω]

allNonvanishingOddMatrixElements[{n1_Integer, m1_Integer}, 
                                 {α_, β_, γ_}, ω_] := 
uniteElements @ nonvanishingOddMatrixElements[{n1, m1}, {α, β, γ}, ω]

Here is a quick check of the just-calculated nonvanishing matrix elements. We compare them with the values of the matrix elements calculated by straightforward integration for explicit values of (Subscript[n, 1]) and (Subscript[m, 1]). To check that we did not miss any nonvanishing matrix elements, we calculate explicitly all matrix elements in the extended convex hull of the nonvanishing {n,m} points.

(* generalized Hénon-Heiles Hamiltonian *)
ℋOp[r_, φ_] = (-ℏ^2/(2 M) (1/r D[r D[#, r], r] + 1/r^2 D[#, φ, φ]) + 
               (α r^4 + β r^2 + γ r^3 Cos[3 φ]) #)&;

checkMatrixElements[{n1_Integer, m1_Integer}, evenOdd: (Even | Odd)] :=
Module[{mF, ψ, m0, mEsS, mEsE},
(* even or odd? *)
Which[evenOdd === Even, 
       mF = allNonvanishingEvenMatrixElements; ψ = ψEven; m0 = 0,
       evenOdd === Odd,  
       mF = allNonvanishingOddMatrixElements; ψ = ψOdd; m0 = 1];
 (* symbolically calculated nonvanishing matrix elements *)
 mEsS = Factor //@ Together //@ mF[{n1, m1}, {α, β, γ}, ω];
 (* explicitly calculated nonvanishing matrix elements *)
 mEsE = Block[{ℏ = 1, M = 1}, Factor //@ Together //@ 
  (Function[{n2, m2}, {n2, m2} -> 
    Integrate[ψ[n1, m1, r, φ] (ℋOp[r, φ] @ ψ[n2, m2, r, φ]) r,
               {φ, 0, 2Pi}, {r, 0, Infinity}, Assumptions -> ω > 0]] @@@
 (Flatten[(* all n and m values covering nonvanishing mEsS values *)
 Table[{n, m}, {n, Max[0 , #1[[1]] - 1], #1[[2]] + 1}, 
               {m, Max[m0, #2[[1]] - 1], #2[[2]] + 1}]& @@
         ({Min[#], Max[#]}& /@ Transpose[First /@ mEsS]), 1]))];     
 (* are nonvanishing elements identical? *)
 Sort[mEsS] === Sort[DeleteCases[mEsE, _ -> 0]]]

checkMatrixElements @@@ 
            {{{0, 3}, Even}, {{1, 1}, Even}, {{2, 2}, Odd}}

After having evaluated the function numberEvenBasisFunctions (or numberOddBasisFunctions respectively), the variable dim has a value and the function nmNumbers has downvalues. The function HillMatrixElements generates a list of all nonvanishing Hill matrix elements for the currently (through nmNumbers) active quantum numbers (Subscript[n, 1]) and (Subscript[m, 1]).

HillMatrixElements[{α_, β_, γ_}, ω_] := 
With[{(* even or odd basis? *)
      anme = If[nmNumbers[1] === {0, 0}, 
                allNonvanishingEvenMatrixElements,
                allNonvanishingOddMatrixElements]},
(* all nonvanishing matrix elements *)
DeleteCases[Flatten[
Table[({k, efNumber[#1]} -> #2)& @@@ anme[nmNumbers[k], {α, β, γ}, ω],
                       {k, dim}]], {_, _efNumber} -> _]];

At this point we also change the origin of the potential minimum. We will set the potential minimum to have the value zero. We do this to make sure that all eigenvalues are strictly positive, so that the function Eigenvalues with a negative integer second argument applied to the resulting sparse array will always yield the lowest eigenvalues.

(* change diagonal matrix elements so that the 
   potential minimum has value zero *)
gaugePotentialMinimum[sa_, {α_, β_, γ_}] :=
Module[{Vmin},
       (* global minimum of the potential *)
       Vmin = Minimize[{α (x^2 + y^2)^2 + β (x^2 + y^2) + 
                        γ (x^3 - x y^2)}, {x, y}][[1]] // 
                                     N[#, Precision[sa]]&;
       (* global gauge transformation to zero potential minimum *)
       sa + SparseArray[Table[{i, i} -> -Vmin, {i, Length[sa]}]]]

SparseHillMatrix[{α_, β_, γ_}, ω_] := 
      gaugePotentialMinimum[#, {α, β, γ}]& @
           SparseArray[HillMatrixElements[{α, β, γ}, ω]]

To visualize the sparse Hill matrix, we implement a function sparseArrayGraphics that displays the nonvanishing matrix values. Small values are shown in red, and large values are shown in blue.

sparseArrayGraphics[sa_, opts___] :=
Module[{(* absolutely scaled matrix elements *)
        ars = (#1 -> ArcTan[Abs[#2]]/(Pi/2))& @@@ 
                         Most[ArrayRules[sa]], min, max},
       (* scale matrix elements relative *)
       {min, max} = {Min[#], Max[#]}& @ (Last /@ ars);
       Graphics[{RGBColor[min/#2, 0, #2/max], 
                 Rectangle[#1 - 1/2, #1 + 1/2]}& @@@ ars, 
       Frame -> True, PlotRange -> All, AspectRatio -> Automatic]]

The next input just reproduces the 2D harmonic oscillator eigenvalues (α=γ=0).

numberEvenBasisFunctions[10];
SA = With[{ω = 1.}, SparseHillMatrix[{0, ω^2/2, 0.}, ω]];
SA // Eigenvalues // Reverse

The next graphic shows the nonvanishing matrix elements for the genuine anharmonic case α=γ=1.

numberEvenBasisFunctions[20];
SA = With[{ω = 1}, SparseHillMatrix[{1, ω^2/2, 1}, ω]];
Show[sparseArrayGraphics[SA]]

By construction, the resulting matrix is symmetric.

SA - Transpose[SA] // Normal // Flatten // Union

We have a quick qualitative look at the eigenvectors at this point. We diagonalize the just-visualized matrix with high precision.

evecs = Eigenvectors[N[SA, 60]];

The following are the angular momentum quantum numbers m of the basis states having a nonvanishing contribution to the resulting eigenstate. We see that in addition to the already above used symmetry, we could further split the states into two subgroups (corresponding to their parity with respect to the rays in the directions {cos(2 π/3), ±sin(2 π/3)}). One subgroup contains all states with m=3k, k∈ℕ, the other the remaining states [✶]. This would allow to use still smaller matrix sizes to obtain the same-quality solution. Doing this is straightforward and will not do it here.

(First /@ Cases[#, {_, 1}])&  /@ 
MapIndexed[{nmNumbers[#2[[2]]][[2]], If[Abs[#1] < 10^-30, 1, 0]}&, 
           Take[#, 24]& /@ Take[evecs, -12], {2}]

The next graphic shows the dependence of the ten lowest eigenvalues on ω. In this concrete example, the eigenvalues are nearly independent of ω over a relatively broad range of ω-values 1/2≲ω≲7/2.

(* prepare matrix *)
numberEvenBasisFunctions[30];
SAω[ω_] = N @ SparseHillMatrix[{1, -2, -3}, ω];
SAωN[ω_] = Normal[SAω[ω]];
(* show eigenvalue dependence on ω *)
Show[Graphics[Line /@ Transpose[(* use various ω-values *)
  Table[{ω, #}& /@ Take[Eigenvalues[SAωN[ω]], -10], 
        {ω, 0.1, 10, 0.1}]]], 
     Frame -> True, PlotRange -> {0, 35}]



We repeat the last calculations and investigate the magnitudes of the eigenvector components. Because the components are not smooth functions of the index, we do not display the components directly. Instead we fit a linear function to the logarithm of selected components (the selected components are the largest components starting with the highest states). This graphics that for ω≈3 we have fastest decay of the eigenfunction components for the lowest states.

Module[{evecs}, 
Show[Graphics[(* color each state differently *)
MapIndexed[{Hue[0.8 (#2[[1]] - 1)/10], PointSize[0.006], Point /@ #}&, 
Transpose[Table[evecs = Abs[Eigenvectors[SAωN[ω]]];
(* fit linear function to logarithm of eigenvector components *)
Table[{ω, Coefficient[Fit[
      (* select maximal components; start from the end *)
      Most[Log[Reverse[FoldList[Max, 0, Reverse[evecs[[-j]]]]]]],
                          {1, x}, x], x]}, {j, 10}], 
      {ω, 0.1, 10, 0.1}]]]], Frame -> True, PlotRange -> All]]

Next, we calculate the Hill matrix for a symbolic value of γ and show the dependence of the lowest eigenvalues as a function of γ. The spectrum of γ is identical to the spectrum for -γ because the sign change of γ in γ cos(3φ)(r^3) effectively just flips that potential around the y-axis. Some even and odd states are degenerate [✶], [✶].

Module[{εMax = 24, δγ = 1/36, SAγ},
Show[Graphics[
(* calculate lowest eigenvalues for various γ-values *)
Function[{eo, col}, eo[εMax];
 Table[(* form matrix and calculate lowest eigenvalues *)
       SAγ = Normal[SparseHillMatrix[{1, 0, γ}, 1]];
       (* even states in red; odd states in blue *)
       {col, PointSize[0.005], Point[{γ, #}]}& /@ 
               Take[Eigenvalues[N @ SAγ], -10], 
       {γ, -3, 3, δγ}]] @@@ 
  (* even and odd states *)
  {{numberEvenBasisFunctions, RGBColor[1, 0, 0]},
   {numberOddBasisFunctions,  RGBColor[0, 0, 1]}},
     PlotRange -> All, Frame -> True]]]

To calculate explicit position representations of the eigenfunctions, we need the numerical values of the Laguerre polynomials. To calculate them quickly (meaning using machine arithmetic), we need a numerically stable recursion relation. To avoid the generation of large values for increasing values of angular quantum numbers, we do not calculate the numerical values of (Subscript[L, n]^λ)(z) directly, but rather of the normalized version (Subscript[χ, n]^λ)(z)=n! /(n+λ)!(Subscript[L, n]^λ)(z). The (Subscript[χ, n]^λ)(z) obey the following recursion relation which is numerically stable when run downwards with respect to λ.

(Subscript[χ, n]^λ)(z)=(z+λ+1) (Subscript[χ, n]^(λ+1))(z)-z (n+λ+2) (Subscript[χ, n]^(λ+2))(z)

Here is a quick check of the quality of this recursion. To avoid large polynomial values for larger r, we also incorporate the exponential term exp(-z) from the 2D harmonic oscillator wave function.

Module[{Lz = 25, δz = 1/2, n = 150, Λ = 150,
        χ, zList, LaguerreLExpList},
 (* scaled LaguerreL functions *)
 χ[n_, λ_, z_] := n!/(n + λ)! LaguerreL[n, λ, z];
 (* z-values *)
 zList = N @ Table[z, {z, 0, Lz, δz}];
 (* starting values of the recursion relation *)
 Do[LaguerreLExpList[n, λ] = N[χ[n, λ, N[zList, 100]] Exp[-zList]], 
    {λ, Λ - 1, Λ}];
 (* carry out recursion downwards *)
 Do[LaguerreLExpList[n, λ] = 
     (1 + zList + λ) LaguerreLExpList[n, 1 + λ] - 
     zList (2 + λ + n) LaguerreLExpList[n, 2 + λ], {λ, Λ - 2, 0, -1}];
 (* display logarithm of relative error as a function of λ *)
 ListPlot[Table[{λ, Log[10, #]}& @ Abs[Max[
      Cases[LaguerreLExpList[n, λ]/
            (χ[n, λ, N[zList, 200]] Exp[-zList]) - 1, _?NumberQ]]],
      {λ, 0, Λ - 2}], PlotRange -> All, Frame -> True, Axes -> False]]

After having evaluated numberEvenBasisFunctions or numberOddBasisFunctions, the function nmNumbers has downvalues. For the {n,m}-values contained in the domain of the function nmNumbers, we will calculate numerical values for discretizations radially and azimuthally. We carry out the radial discretization in [0,rMax] and the azimuthal discretization in [0,π] (because we will later visualize the squares of the eigenfunctions and all eigenfunctions have definite parity with respect to the x-axis, we do not need the [π,2π] domain).

cacheEigenFunctionLists[ω_, rMax_, {ppr_, ppφ_}] :=
Block[{χ, n, Λ, LaguerreLExpList, δr = rMax/ppr, L = 1/Sqrt[ω],
       Λming, Λg, tpd = Developer`ToPackedArray},
 Clear[RList, ΦList];
 (* extremal n, m values *)
 nAndMinmMaxmList = {#1[[1, 1]], Min[Last /@ #], Max[Last /@ #]}& /@ 
       Split[Sort[nmNumbers /@ Range[dim]], #1[[1]] == #2[[1]]&];
 Do[usednmPair[nmNumbers[k]] = True, {k, dim}];
 (* radial eigenfunctions *)
 (* list of scaled, squared radii *)
 zList = N @ Table[(r/L)^2, {r, 0, rMax, δr}];
 χ[n_, λ_, z_] := n!/(n + λ)! LaguerreL[n, λ, z];
 Do[{n, Λmin, Λ} = nAndMinmMaxmList[[j]]; 
    (* calculate normalized Laguerre polynomials *)
    Do[LaguerreLExpList[n, λ] = χ[n, λ, zList] Exp[-zList/2], 
       {λ, Λ - 1, Λ}];
    Do[LaguerreLExpList[n, λ] = 
       (1 + zList + λ) LaguerreLExpList[n, 1 + λ] - 
       zList (2 + λ + n) LaguerreLExpList[n, 2 + λ], 
     {λ, Λ - 2, Λmin, -1}];
   (* cache the radial wavefunction values *)  
   Do[If[usednmPair[{n, λ}], 
         (* list of radial wavefunction values *)
         RList[n, λ] = tpd[Sqrt[2]/L 1/Sqrt[n!/(n + λ)!]*
                           If[λ =!= 0, zList^(λ/2),
                              Table[1., {r, 0, rMax, δr}]]*
                               LaguerreLExpList[n, λ]]],
      {λ, Λmin, Λ}], {j, Length[nAndMinmMaxmList]}];
 (* azimuthal eigenfunctions in [0, Pi] *)
 φList = N @ Table[φ, {φ, 0, Pi, Pi/ppφ}]; 
 (* all occurring m quantum numbers *)
 {Λming, Λg} = {Min[#], Max[#]}& @ Rest[Transpose[nAndMinmMaxmList]];
 φEF = If[Λming == 0, ΦEven, ΦOdd];
 Do[ΦList[m] = tpd[N @ φEF[m, φList]], {m, Λming, Λg}]]

The function cacheEigenFunctionLists uses machine precision and will work well with (Subscript[n, max]), (Subscript[m, max]) a few hundred. While this will be large enough for the visualizations to be presented at the later, for very high states, one might have to use high-precision arithmetic to calculate the radial eigenfunctions (the Hill eigenvectors are good enough with machine arithmetic).

The following graphic shows a comparison between the directly calculated radial eigenfunctions R[8, 8, r] (thick gray line) and the discretized eigenfunctions (thin red line) obtained through the function cacheEigenFunctionLists.

Show @
Block[{n = 24, m = 24, ℏ = 1, M = 1, ω = 2, 
       $DisplayFunction = Identity, ppr = 100, rMax = 10},
numberEvenBasisFunctions[100];
cacheEigenFunctionLists[ω, rMax, {ppr, 1}];
{(* exact radial eigenfunctions *)
 Plot[Evaluate[R[n, m, r]], {r, 0, rMax}, PlotPoints -> 200,
      PlotStyle -> {GrayLevel[0.5], Thickness[0.02]}],
 (* recursively-calculated radial eigenfunctions  *)
 ListPlot[MapIndexed[{(#2[[1]] - 1)/ppr rMax, #1}&, RList[n, m]], 
          PlotRange -> All, PlotJoined -> True, 
          PlotStyle -> {Hue[0]}]}]

We visualize the squares of the eigenfunctions using 3D plots and contour plots. Because the discretized eigenfunctions are defined on a tensor product grid in polar coordinates, we implement the two functions polarListPlot3D and polarListContourPlot that display such type discretized functions over a circular domain. The argument of both function is the matrix gef of function values in polar coordinates in the upper half-plane.

polarListPlot3D[gef_, rMax_:1, color___] :=
Module[{ppr, ppφ, r, φ, points, polysU, polysL},
 {ppr, ppφ} = Dimensions[gef] - 1;
 (* Cartesian coordinates of points *)
 points = MapIndexed[({r, φ} = {rMax (#2[[1]] - 1)/ppr, 
                                Pi (#2[[2]] - 1)/ppφ};
                      {r Cos[φ], r Sin[φ], #1})&, gef, {2}] // N;
 (* polygons in upper half-plane *)
 polysU = Table[Polygon[{points[[i, j]], points[[i + 1, j]],
                         points[[i + 1, j + 1]], points[[i, j + 1]]}], 
                {i, ppr}, {j, ppφ}];
 (* polygons in lower half-plane *)
 polysL = Map[(# {1, -1, 1})&, polysU, {-2}];              
 (* resulting 3D graphics *)
 Show @ Graphics3D[{EdgeForm[], color, polysU, polysL}, 
      PlotRange -> All, Boxed -> False, BoxRatios -> {1, 1, 0.4}]]

To see the structure in the squares of the higher-order wave functions better, we use a nonlinear scaling for the coloring of the contour regions. Actually, we implement two color functions. One is based on a nonlinear function applied directly to the function values and the second on a linear function applied on the number of function values. The first one shows the overall structure of the eigenfunctions (having physical relevance), while the second reveals fine details (showing numerical details, largely resulting from the summation) in the regions where the function values are small.

polarListContourPlot[gef_, resolveFineStructure_: False, rMax_:1] :=
Module[{ppr, ppφ, λ, 𝒸𝒻, cts = 48},
 {ppr, ppφ} = Dimensions[gef];
 If[resolveFineStructure,
    (* homogeneously spaced contours *)
    cls = Union[#[[Round[ppr ppφ/cts/2]]]& /@ 
                Partition[Sort[Flatten[gef]], Round[ppr ppφ/cts]]];
    λ = Length[cls];
    (* use linear gradient of graylevels with respect to contours *)
    𝒸𝒻 = Function @@ {Which @@
      Flatten[Append[MapIndexed[{#1, GrayLevel[1 - (#2[[1]] - 1)/λ]}&, 
                     Function[v, # < v] /@ cls], {True, GrayLevel[0]}]]};
    colorOpts = Sequence[ColorFunctionScaling -> False,
                         ColorFunction -> 𝒸𝒻, Contours -> cls],
    (* nonlinear color function scaling *)
    𝒸𝒻[x_] :=  GrayLevel[1 - If[x < 1/2, 0 + Sqrt[1/4 - (x - 1/2)^2], 
                                         1 - Sqrt[1/4 - (x - 1/2)^2]]];    
    colorOpts = Sequence[ColorFunctionScaling -> True,
                         ColorFunction -> 𝒸𝒻, Contours -> cts]];                  
 (* contour plot in the upper half-plane *)
 lcp = ListContourPlot[gef, colorOpts, 
            DisplayFunction -> Identity, ContourLines -> False,
            PlotRange -> All, MeshRange -> {{0, Pi}, {0, rMax}}, 
            ContourStyle -> {{GrayLevel[0], Thickness[0.002]}}];
 (* make a graphics and map polygons to Cartesian coordinates *)
 gr = Graphics[lcp];
 grR = gr /. (pl:(Polygon | Line))[points_] :>   
              pl[addPoints[points, 10^-2 Min[rMax, 2Pi]]];
 (* upper half-plane polygons *)
 grRCU = grR /. (pl:(Polygon | Line))[points_] :>   
                 pl[{#2 Cos[#1], #2 Sin[#1]}& @@@ points];  
 (* lower half-plane polygons *)
 grRCL = grRCU /. (pl:(Polygon | Line))[points_] :>   
                   pl[(# {1, -1})& /@ points];                
 (* resulting contour graphic *)                
 Show[{grRCU, grRCL},
      Frame -> False, DisplayFunction -> $DisplayFunction]]

To generate a contour plot in polar coordinates, we map a contour plot made in "Cartesian" φ,r-coordinates to Cartesian x,y-graphics coordinates. Because the polygons used to represent a contour plot have partially long edges, we have to subdivide these edges to obtain a smooth polygon boundary after the mapping to Cartesian graphics coordinates. The function addPoints inserts points into the list of points points, so that the largest resulting line segment is smaller than δε.

(* add intermediate points into list points so that the maximal
   distance between two consecutive points is δε *)
addPoints[points_, δε_] := 
 Module[{n, l}, Join @@ 
       (Function[s, (* segment s short enough? *)
               If[(l = Sqrt[#. #]&[Subtract @@ s]) < δε, s, 
         n = Floor[l/δε] + 1; (* form segments *)
         Table[# + i/n (#2 - #1), {i, 0, n - 1}]& @@ s]] /@ 
              (* form segments between points *)
              Partition[Append[points, First[points]], 2, 1])]

To construct the position representation of the eigenfunctions, we use the function constructEigenFunction. Its argument is an eigenvector from the diagonalized Hill matrix. The function constructEigenFunction assumes that the functions numberEvenBasisFunctions and cacheEigenFunctionLists have already been evaluated, so that the current basis functions are globally visible. We form the tensor product of the azimuthal and the radial eigenfunctions for only those states, contributing significantly to the sum (here we use the value (10^(-4)) for the relative contribution, which is sufficient for graphics purposes).

constructEigenFunction[evec_] :=
Module[{maxBasisValues, relWeightList, εR, ef, n, m,
        (* small enough for graphics purposes *) ε = 10^-4},
       (* relative size of products of expansion coefficients 
          times [radial part of the] eigenvector magnitude *)
       maxBasisValues = Table[Max[Abs[RList @@ nmNumbers[j]]], 
                              {j, Length[evec]}];
       relWeightList = maxBasisValues Abs[evec];
       εR = ε Max[relWeightList];                     
       (* assume RList is a global variable *)
       ef = Table[0., {Length[RList[0, 1]]}, {Length[ΦList[1]]}];
       (* add all eigenfunctions that contribute significantly *)
       Do[If[relWeightList[[j]] > εR,
             {n, m} = nmNumbers[j]; 
             (* faster method for: 
                ef = ef + evec[[j]] Outer[Times, RList[n, m], ΦList[m]]; *)
             LinearAlgebra`BLAS`GER[Evaluate[evec[[j]]], 
                                       RList[n, m], ΦList[m], ef]],
          {j, Length[evec]}];
       (* return summed eigenfunction *) ef]

We start with a calculation of the states for the parameter values α=1, β=-1/2, and γ=-3. We use a relatively small basis of 650 states (up to energy 50) and diagonalize the full (dense) matrix). Here are the lowest eigenvalues.

(* number all basis function up to energy 50 (basis size 650) *)
numberEvenBasisFunctions[50]; Print["dim = ", dim]; 
(* form sparse Hill matrix *)
SA = With[{ω = 1}, SparseHillMatrix[N @ {1, -1/2, -3}, N @ ω]]; 
(* calculate eigenvalues and eigenvectors *)
esystem = Eigensystem[N[Normal[SA]]]; 
(* sort by increasing eigenvalues *)
aux = Sort[MapIndexed[{#1, #2[[1]]} & , esystem[[1]]]]; 
evecs = esystem[[2, Last /@ aux]];
First /@ Take[aux, 6]

A quick look at the eigenvectors shows that this basis is a bit too small for a reliable visualization (we will use larger bases in a moment). But for a qualitative correct graphics we have enough information.

Show[GraphicsArray[
Table[ListPlot[Log[10, Abs[evecs[[j]]]], 
               PlotRange -> {0, -6}, AxesOrigin -> {0, 0},
               DisplayFunction -> Identity], {j, 3}]]]

Now, we calculate the discretized radial eigenfunctions and display the squares of the eigenfunctions as 3D plots. The qualitative shape of the eigenfunctions is as expected: in the ground state we see the full symmetry of the potential. The next higher states still have the reflection symmetry with respect to the x-axis. The second even state has only a slightly higher energy than the ground state (for three independent wells, the splitting will be exponentially small).

(* calculate position representation of eigenfunctions *)
cacheEigenFunctionLists[1, 6, {60, 120}]

Show[GraphicsArray[
 Block[{$DisplayFunction = Identity}, 
  Table[polarListPlot3D[constructEigenFunction[evecs[[j]]]^2, 6,
                        SurfaceColor[Hue[0.12], Hue[0], 2.5]], 
        {j, 3}]]]]

Here is a higher state shown.

polarListPlot3D[constructEigenFunction[evecs[[24]]]^2, 6,
                        SurfaceColor[Hue[0.78], Hue[0.12], 2.26]]

Here are contour plots of some of the higher states.

Show[GraphicsArray[
 Block[{$DisplayFunction = Identity}, 
  Table[polarListContourPlot[constructEigenFunction[evecs[[j]]]^2, False], 
        {j, 10, 12}]]]]

Using True for the second argument of polarListContourPlot yields plots whose fine structure reflects the details of the numerical evaluation of the eigenfunctions. Changing the basis size of the azimuthal and radial discretization would show differently-looking plots for the setting True.

(DownValues[In, Sort -> True][[-2]] /. False -> True)[[2]]

Now, we visualize some of higher states for some selected potentials. Due to their considerably more complicated structure, the contour plot generated by polarListContourPlot is frequently easier to grasp than the corresponding 3D plot from polarListPlot3D. Given a parity evenOdd, set of potential parameters β and γ, the maximal energy of the 2D harmonic oscillator basis states ε, and its frequency ω, the function eigenFunctionContourPlotArray calculates and displays ℴ states around the energy shift above the potential bottom. The option showEigenvectorComponents shows the base 10 logarithm of the absolute values of the eigenvector component. This allows to judge the correctness of the contour plots.

(* estimate spatial extension of a state with energy ε 
   in a potential with parameters {α, β, γ} *)
rClassical[ε_, {α_, β_, γ_}] := Max[Function[pm1, 
Select[ρ /. {ToRules[NRoots[ε == α ρ^4 + β ρ^2 + pm1 γ ρ^3, ρ]]},
       (Im[#] == 0 && # > 0)&]] /@ {+1, -1}]

eigenFunctionContourPlotArray[evenOdd: (Even | Odd), 
                              {β_, γ_}, {εMax_, ω_}, {shift_, ℴ_},
                              {ppr_, ppφ_}, opts___] :=
Module[{SA, evecs, ℛ, gef, λ = 3}, 
Print[{β, γ}];
If[evenOdd === Even, numberEvenBasisFunctions, 
                     numberOddBasisFunctions][εMax]; 
Print["dimFinal = ", dim]; 
(* use high precision to guarantee that matrix is symmetric *)
SA = N[SparseHillMatrix[{1, β, γ}, SetPrecision[ω, 30]]]; Print[SA];
(* calculate eigenvalues *)
evecs = Reverse @ 
Eigenvectors[SA, ℴ, Method -> {"Arnoldi", Tolerance -> 10^-8,
             BasisSize -> Min[dim, Abs[20 ℴ]], MaxIterations -> 10^4,
             Shift -> shift}]; Print["evecs done"];
(* approximate spatial extension of the wavefunction *)            
ℛ = 1.1 rClassical[N[shift], {1, β, γ}] + 1/4; Print["ℛ = ", ℛ];
(* calculate discretized 2D HO eigenfunctions *)  
cacheEigenFunctionLists[ω, ℛ, {ppr, ppφ}];
(* partition a list into length λ; append potential rest *)  
partition[l_] := If[Mod[Length[l], λ] === 0, Partition[l, λ],
         Append[Partition[l, λ], Take[l, -Mod[Length[l], λ]]]];
(* potentially show eigenvectors *)  
If[(# /. {opts} /. # -> False)&[showEigenvectorComponents],
   Show[GraphicsArray[#]]& /@ 
   partition[ListPlot[Log[10, Abs[#]], PlotRange -> {-6, 0}, 
                      DisplayFunction -> Identity, 
                      AxesOrigin -> {0, 0}]& /@ evecs]];
(* show eigenfunctions squared *)     
Show[GraphicsArray[#]]& /@ partition[
Table[gef = constructEigenFunction[evecs[[j]]];
      Block[{$DisplayFunction = Identity}, 
            polarListContourPlot[Abs[gef]^2]], {j, ℴ}]]]

Now we have all functions implement to lean back and enjoy a gallery of Hénon–Heiles potential eigenfunctions. While the function eigenFunctionContourPlotArray displays list plots of the Hill matrix eigenvectors, we do not show them explicitly. The reader can readily rerun the code to check that all eigenvectors decay rapidly enough to allow for a reliable visualization. We will always display six energetically consecutive states. In average, every third state has the full threefold rotational symmetry.

We start with some relatively simple eigenfunctions using 6480 basis functions. The reflection symmetry on the real axis and the relicts of the threefold rotational symmetry are clearly visible. We also show the size of the eigenvectors of the Hill matrix. The eigenvectors decay quickly. The decay of the Hill matrix eigenfunctions depends on the chosen ω. As a rule of thumb, we use in the following a value of ω, such that (Subscript[r, c])≈k ω where (Subscript[r, c]) is the maximal classical extension of a state and k is a small integer.

eigenFunctionContourPlotArray[Even, {205/42, -13/3}, {160, 3/2}, {57, 6},
               {90, 120}, showEigenvectorComponents -> False]

Here are some slightly more complicated states. We use 8190 basis functions. Neighbor states can have quite different-looking probability distributions.

eigenFunctionContourPlotArray[Even, {-49/47, 149/86}, {180, 3/2}, {106, 6},
               {90, 120}, showEigenvectorComponents -> False]

And the next set of states shows only a minor influence of the three wells on the probability distributions.

eigenFunctionContourPlotArray[Even, {487/100, 191/31}, {220, 3/2}, {101, 6},
                  {120, 160},  showEigenvectorComponents -> False]

We continue with two sets of higher states having more structures.

eigenFunctionContourPlotArray[Even, {70/17, -249/55}, {200, 3/2}, {95, 6},
                          {90, 120}, showEigenvectorComponents -> False]

For the next states, we use more plot points to obtain a better resolved graphics of their fine structure.

eigenFunctionContourPlotArray[Even, {259/44, 117/47}, {180, 2}, {180, 6},
                    {120, 160}, showEigenvectorComponents -> False]

To visualize still higher states reliably, we use a larger basis set containing 14520 functions.

eigenFunctionContourPlotArray[Even, {247/88, 163/34}, {240, 2}, {210, 6},
                        {120, 160}, showEigenvectorComponents -> False]

We now move on to the odd states. The odd states look qualitatively similar to the even states. Here are again some relatively simple states.

eigenFunctionContourPlotArray[Odd, {-481/85, -17/8}, {160, 3/2}, {95, 6},
                   {90, 120}, showEigenvectorComponents -> False]

Here are three sets of slightly more complicates states.

eigenFunctionContourPlotArray[Odd, {-5, -19/6}, {180, 2}, {145, 6},
                {120, 160}, showEigenvectorComponents -> False]

eigenFunctionContourPlotArray[Odd, {286/83, 8/3}, {180, 2}, {155, 6},
                      {120, 160}, showEigenvectorComponents -> False]

eigenFunctionContourPlotArray[Odd, {479/80, 423/97}, {180, 2}, {133, 6},
                   {120, 160}, showEigenvectorComponents -> False]

We continue with some more complicated states. We again use a larger basis size.

eigenFunctionContourPlotArray[Odd, {-8/11, -36/11}, {180, 5/2}, {171, 6},
                      {120, 160}, showEigenvectorComponents -> False]

We end with some still higher states of the potential from the last example.

eigenFunctionContourPlotArray[Odd, {-8/11, -36/11}, {200, 2}, {220, 6},
                      {150, 200}, showEigenvectorComponents -> False]