The EKB quantization

Quantum Mechanics
                      K. Konishi, G. Paffuti

Study the EKB quantization procedure for the potentials

V[x,y]= (1/2)((Subscript[ω, x])*(x^2)+(Subscript[ω, y])*(y^2) )+ λ x((y^2)+ η (x^2));

V[x,y]= ((ω^2)/2)((x^2)+(y^2) )+ (g/2) (x^2)((y^2)+  (x^2)).

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Introduction

In the complements to the chapter on WKB approximation it has been shown that the extension to higher dimensions, for integrable systems, of the semiclassical quantization procedure is given by the EKB condition

(1/(2π ℏ))(∮_(Subscript[γ, i])) (∑_s) (Subscript[p, s])d(Subscript[q, s]) = (n+∑(α/4))

where γ is a non contractible path on the invariant torus defined by integral of motions, and the sum runs over the crossings of γ with caustics. For usual reflection points α = 1.

The previous formula can be extended to small perturbation of an integrable system, in situations where the existence of invariant tori is granted by the KAM theorem.

For two dimensional systems a possible strategy to implement practically this recipe is to use Poincare sections: take the projection of a trajectory on the x - (Subscript[p, x]) plane. If,  as it is usual, the trajectory fill invariant torus, his projection traces a closed curve in this plane. In this plane dy = 0 then the quantization condition reduces to

(1/(2π ℏ))(∮_γ)  (Subscript[p, x])dx = ((Subscript[n, x])+(1/2))

The same procedure can be used for the projection on y-(Subscript[p, y]) plane, with a corresponding quantum number (Subscript[n, y]). Varying energy and a second parameter we can select trajectories which satisfy these constraints, obtaining a semi-classical approximation for the eigenvalues.

In the following we take as parameters to explore torus quantization the energy and a parameter α connected to (Subscript[p, x ])initial momentum with

(Subscript[p, x])[0] = ((Sqrt[2 K])/(Sqrt[1+(α^2)]))

K is the initial kinetic energy. We use α because can be freely varied, independently on energy.

We present to applications:

Functions

Let us suppose we have defined a potential V[x,y] and an initial point (x0,y0) for our trajectories. The function getPoincareSections below compute Poincaré sections, giving as outputs two lists with the couples (x,(Subscript[p, x]))and (y,(Subscript[p, y])) respectively. The first are x and relative momentum of the trajectory in the intersection with y=0 plane. The second list is the analogous couple for intersections with x = 0 plane.

The technique for getting Poincaré sections is exposed in the Mathematica  help and use the EventLocator method in NDSolve. The Event in our function is given by the couple {y[t], x[t]} i.e. the function detect times for which y[t]=0 or x[t] = 0. The option Direction-> {1,1} select events for growing values of the coordinates, i.e. zeros for passages from negative to positive values.

The EventAction uses the Reap-Sow mechanism of Mathematica. We select couples {x[t],px[t]} or {y[t],py[t]} and assign to them to different tags, tagx, tagy. The function Reap gives as output these two lists, separated with the tags, passed as arguments of Reap. The last element of the output of Reap contains the two lists. These are selected with

{dxpx, dypy} = outReap[[-1]]outReap[[-1]].

It is essential to use a method which preserve energy in the numerical solution. We use SymplecticPartitionedRungeKutta, a submethod of NDSolve (see Mathematica Help).

getPoincareSections[ene_,α_,tmax_]:=
 Module[{eneKin,x,y,px,py,px0,py0,t,eqs,ics,dxpx,dypy,outReap,tagx,tagy},
eneKin = ene - V[x0,y0];
If[eneKin ≤ 0,Print["energy too small"];Abort[];];
px0 = ((Sqrt[2 eneKin])/(Sqrt[1+(α^2)])); py0 = (Sqrt[2 eneKin - (px0^2)]);
sprkmethod={"SymplecticPartitionedRungeKutta", "DifferenceOrder"->4,"PositionVariables"-> {x[t],y[t]}};
eqs= {x'[t] ==px[t],y'[t]== py[t],px'[t]== - D[V[x[t],y[t]],x[t]],
py'[t]== - D[V[x[t],y[t]],y[t]]};
ics= {y[0]==y0,py[0]==py0,x[0]==x0,px[0]==px0};
outReap=Reap[NDSolve[{eqs,ics},{},{t,0,tmax},
Method->{"EventLocator","Event":>{y[t],x[t]},
"Direction"->{1, 1},"EventAction":>If[Abs[y[t]]≤ (10^(-6)),Sow[{x[t],px[t]},tagx],Sow[{y[t],py[t]},tagy]],Method->sprkmethod},
MaxSteps->∞],{tagx,tagy}];
{dxpx,dypy}=outReap[[-1]];
{dxpx[[1]],dypy[[1]]}]

Once a Poincaré section is obtained we have to compute the action integral for this curve. We use a simple method which works if the curve is convex, it must be changed in the case of non convex curves. We use the symmetry p → - p and select points with positive and negative momenta. if the curve is convex each set contains a list (q,p) which can be used to construct an interpolatng function which can be integrated with Integrate. The function integralP compute the integrals, while action[list_] accepts the lists (x,(Subscript[p, x])) or (y,(Subscript[p, y])) and after a couple of calls to integralP make the difference of the two integrals to construct the integral on the closed curve.

action[listXP_]:= Module[{pPos,pNeg},
{pPos,pNeg}={Select[listXP,#1[[2]]≥ 0&],Select[listXP,#1[[2]]< 0&]};
 (1/(2π))(integralP[pPos]-integralP[pNeg])]

actionProva[listXPIn_]:= Module[{pPos,pNeg,listXP,mean},
If[switchXP == False, mean = 0;
listXP=listXPIn;,
listXP= Transpose[{listXPIn[[All,2]],listXPIn[[All,2]]}];
 mean = Mean[listXP[[All,2]]];
];
{pPos,pNeg}=
{Select[listXP,#1[[2]]≥ mean&],Select[listXP,#1[[2]]< mean&]};
(* (1/(2π))(integralP1[pPos]+integralP1[pNeg]) *)
 (1/(2π))(integralP[pPos]-integralP[pNeg])]

integralP[listXP_]:= Module[{list1,xmin,xmax,xpOrdered,pOfx},
xmin = Min[listXP[[All,1]]];
xmax = Max[listXP[[All,1]]];
xpOrdered= Sort[listXP,#1[[1]]≤ #2[[1]]&];
pOfx = Interpolation[xpOrdered,InterpolationOrder-> 2];
Integrate[pOfx[x],{x,xmin,xmax}]]

The function valuesNxNy just compute the two actions in x and y variables in a form compatible with FindRoot. FindRoot is used in the function quantize to impose quantization conditions. We take as an initial guess for eigenvalues those of the uncoupled modes, i.e. the second derivative of the potential in (0,0) (which is supposed to be a local minimum).

valuesNxNy[ene_?NumberQ,α_?NumberQ]:= Module[{lis,nx,ny,dxpx,dypy},
{dxpx,dypy}= getPoincareSections[ene,α,tmax];
nx=action[dxpx];
ny=action[dypy];
{nx,ny}];

quantize[nx_,ny_]:= Module[{ω1,ω2,eneIn,en,α,x,y},
ω1 =Sqrt[( (∂_(x,x))V[x,y]/.{x-> 0,y-> 0})];
ω2 =Sqrt[( (∂_(y,y))V[x,y]/.{x-> 0,y-> 0})];
eneIn = (nx+0.5)ω1+(ny+0.5)ω2;
{en,α}/.FindRoot[valuesNxNy[en,α]== {nx+0.5,ny+0.5},{{en,0.9eneIn,1.1eneIn},
{α,0.5,1.5}},PrecisionGoal-> 4,
MaxIterations-> 40]]

Finally the function getTrajectory compute the trajectory for fixed values of the parameters, this can be useful for plotting purposes.

getTrajectory[ene_,α_,tmax_]:= Module[
{eneKin,x,y,px,py,px0,py0,t,eqs,ics,dxpx,dypy,pap,tagx,tagy},
eneKin = ene - V[x0,y0];
If[eneKin ≤ 0,Print["energy too small"];Abort[];];
px0 = ((Sqrt[2 eneKin])/(Sqrt[1+(α^2)])); py0 = (Sqrt[2 eneKin - (px0^2)]);
sprkmethod={"SymplecticPartitionedRungeKutta", "DifferenceOrder"->4,"PositionVariables"-> {x[t],y[t]}};
eqs= {x'[t] ==px[t],y'[t]== py[t],px'[t]== - D[V[x[t],y[t]],x[t]],
py'[t]== - D[V[x[t],y[t]],y[t]]};
ics= {y[0]==y0,py[0]==py0,x[0]==x0,px[0]==px0};
{x,y,px,py}/.First[NDSolve[{eqs,ics},{x,y,px,py},{t,0,tmax},Method->sprkmethod,MaxSteps->∞]]];

The function hermiteMatrices is reported from problem [*] and it will be used to compute eigenvalues for a quartic oscillator in the variational approach.

hermiteMatrices[dim_,parity_Integer:-1]:=
Module[{x2,p2,x4,x1}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[parity == 0 || parity == 1,Return[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

Cubic Potential

We take as first example the potential

V[x_,y_]:= (1/2)(ωx2*(x^2)+ωy2*(y^2) )+ λ x((y^2)+ η (x^2));

See ref[*]

{ωx2,ωy2,λ,η} = {0.29385,2.12581,-0.1116,0.08414};
x0 = 0; y0=0;

This potential has some metastable states, from ref[*] we get eigenvalues eneTH[ (Subscript[n, x]), (Subscript[n, y])]:

eneTH[0,0] = 0.9916; eneTH[1,0]=1.5164; eneTH[2,0]=2.0308; eneTH[0,1]=2.4188;

We choose the time used for evolution of trajectories

tmax = 1000;

and obtain semiclassical eigenvalues with quantize.

N.B.  The procedure can be quite time consuming as we used no tricks to select parameters. The reader can use, if needed, a smaller value for tmax, or change the PrecisionGoal in the function FindRoot used in quantize.

Off[NDSolve::"noout"];

({ene,frac}=quantize[0,0])//Timing

The result is in good agreement with quantum mechanical numerical result. Let us look at the Poincaré sections:

{dxpx,dypy}=getPoincareSections[ene,frac,tmax];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]

With longer times the sections tend towards continuous lines. It can happens that for particular values of energy the orbit is periodic, in this case the number of points do not change with time and a different procedure has to be used for integrals.

Let us check also other states

{ene,frac}=quantize[1,0]

{dxpx,dypy}=getPoincareSections[ene,frac,2000];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]

{ene,frac}=quantize[0,1]

{dxpx,dypy}=getPoincareSections[ene,frac,2000];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]

All the results are quite close to theoretical values :

eneTH@@#1&/@{{0,0},{1,0},{2,0},{0,1}}

Trajectory

Here we show a plot of the trajectory in ordinary space , the large ellipsoid is the region allowed by energy conservation. Caustics are apparent.

{x,y,px,py}=getTrajectory[ene,frac,2000];

Show[{ParametricPlot[{x[t],y[t]},{t,0,2000},PlotRange-> All],
ContourPlot[V[x1,y1]== ene,{x1,-4,5},{y1,-3,3}]},PlotRange-> All]

It is impossible to draw a threedimensional plot of a four - dimensional torus, then we limit to show the 3D view of the orbit in the variables (x,y,(Subscript[p, x])):

ParametricPlot3D[{x[t],y[t],px[t]},{t,0,2000},PlotRange-> All,
AxesLabel-> {"x","y","(Subscript[p, x])"}]

Quartic potential

As a second example we consider a quartic potential.

V[x_,y_]:= ((ω^2)/2)((x^2)+(y^2) )+ (g/2) (x^2)((y^2)+  (x^2));
x0 = 0.0;y0=0.0;

ω = 1.5; g = 0.15;

Variational computation

First let us compute eigenvalues using the techniques developed in problem [*]. We construct a matrix Hamiltonian using tensor products for matrices in x and y spaces.

dim = 30;
{p2,x2,x4,x1}=N/@hermiteMatrices[dim];

unitM = IdentityMatrix[dim];
a_⊗b_ := KroneckerProduct[a,b]

ham =(1/2) p2⊗unitM + (1/2)unitM⊗p2 +(ω^2)(1/2)x2⊗unitM+(ω^2)(1/2)unitM⊗x2 + (1/2)g(x2⊗x2+ x4⊗unitM);

Let us compute the eigenvalues:

eneVariational=Sort[Eigenvalues[Normal[ham]]][[1;;30]]

EKB quantization

We use the same procedure adopted in the cubic case.

({ene,frac}=quantize[0,0])//Timing

The result agree with variational result within 1%.

These are the Poincaré sections:

{dxpx,dypy}=getPoincareSections[ene,frac,tmax];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]

and this is the corresponding trajectory

{x, y, px, py} = getTrajectory[ene, frac, 1000]; 
ParametricPlot[{x[t],y[t]},{t,0,800}]

Similar results canbe found for first excited states

{ene,frac}=quantize[1,0]

{dxpx,dypy}=getPoincareSections[ene,frac,1000];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]

or for relatively high excited states

{ene,frac}=quantize[5,1]

The coresponding numerical eigenvalue was 11.9036.

{dxpx,dypy}=getPoincareSections[ene,frac,1000];
GraphicsRow[{ListPlot[dxpx,AxesLabel-> {"x","(Subscript[p, x])"}],ListPlot[dypy,
AxesLabel-> {"y","(Subscript[p, y])"}]},ImageSize-> 72 6]