WKB approximation

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study numerically the energy levels and wave functions for harmonic and anharmonic potentials in the WKB approximation.

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

Introduction

We give some simple applications of the WKB methods. The attention is mainly on exact resolvable cases in order to check the methods. For other cases a variational computation is taken as reference, the functions reported in section Functions are those studied in the chapter on variational calculations.

Functions

We report here some functions useful to perform variational computations. For detailed informations on their use see notebooks in the chapter on variational techniques.

The main function is hermiteMatrices which compute in matrix form operators (p^2), (x^2), etc. The calling sequence is

{p2, x2, x4} = hermiteMatrices[dim, parity]

It produces matrix elements for operators p^2 , x^2, x^4. Parity is 0 or 1 for even and odd states. dim is the dimension of the basis. If parity has any other value full matrices are computed, for both even and odd states, and the output is changed adding matrix elements of x in the last place. The default is parity = -1, i.e. full matrices

{p2, x2, x4,x1} = hermiteMatrices[dim, -1],  or  hermiteMatrices[dim].

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[Normal/@{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
Normal/@{p2,x2,x4,x1}]];

Basis' s functions are given in general by

fBase[k_,x_]:= (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,x]Exp[-x^2/2]

For fixed parity:

fBase[n_,x_,par_]:=(1/(π^(1/4)))(1/(Sqrt[(2^(2n-2+par))(2n-2+par)!]))HermiteH[2n-2+par,x]Exp[-x^2/2]

Wave function computation

Here we just write recursion equation for Hermite polynomials, the reader usually will not need a direct call to this function. It is compiled for speed and the integer n represents the number of polynomials, i.e. all polynomials from degree 0 to degree n-1.

recursiveHermitePolynomials= Compile[{{n,_Integer},{x,_Real,1}},
Module[{y1,y0,p1,ris,ris0},
y0 = 1.0 + 0.*x;
y1 = 2.x;
ris0={y0,y1};
ris=Table[
p1=2. x y1 - 2.  k y0;
y0=y1;
y1=p1;
p1,{k,1,n-2}];
ris=Join[ris0,ris];
ris]];

This function computes normalized Hermite polynomials, in the form needed for oscillator wave functions. We write two versions of it, one with arguments [dim, r] gives all polynomials up to degree dim-1.  The second with three arguments, [dim, r, par], gives the first dim polynomials of parity fixed by par, (0 for even polynomials, 1 for odd polynomials).

hermiteP[dim_,xIn_]:= Module[{cf,x,ris,coeff},
x = Flatten[{xIn}];
 coeff=Table[(1/((N[π])^(1/4))) (1/(Sqrt[(2.^(k-1))(k-1)!])),{k,1,dim}]; 
ris= recursiveHermitePolynomials[dim,x]*coeff;
If[Length[xIn]== 0,Flatten[ris],ris]];

hermiteP[dim_,xIn_,par_]:= Module[{cf,dimUp,index,x,ris},
x = Flatten[{xIn}];
dimUp=2 dim - 2 + par+1;
index=Range[par+1,2 dim-2+par+1,2];
ris=hermiteP[dimUp,x][[index]];
ris];

A more efficient version would use recurrence only for even or odd polynomials, we leave to the reader the task of writing such a routine.

Finally getPsi give the wave function for a coefficient list eigv at points r. Parity can be specified, if not all polynomials are considered. This is the default.

getPsi[eigv_,r_,par_: -1]:=
If[par== 1|| par==0, 
(eigv.hermiteP[Length[eigv],r,par])*Exp[-(r^2)/2],
(eigv.hermiteP[Length[eigv],r])*Exp[-(r^2)/2]];

The harmonic oscillator

We know analytically both the exact and the WKB solution for this problem. Using natural units, ℏ = m = 1, the model is defined by a potential

V[x_]:= (1/2)(x^2)

The exact solutions are (n= 0,1,... )

ψThHO[n_,x_]:=  (1/(π^(1/4))) (1/(Sqrt[(2^n)(n)!]))HermiteH[n,x]Exp[-x^2/2]

while the WKB approximated wave functions are (see text and problems in this chapter):

ψHO[n_,x_]:=Module[{a,C},  (* a is the turning point (Sqrt[((2 E)/(m (ω^2)))]) *)
a = (Sqrt[2. (n+1/2)]);
(* C is the normalization constant (Sqrt[((2m ω)/π)]) *)
C = (Sqrt[2/π]);
 Which[
x<0,(((-1))^n)ψHO[n,-x],
x< a, (C/((((a^2)-(x^2)))^(1/4)))Cos[((π n)/2)+(n+(1/2))((x/a) (Sqrt[1-((x^2)/(a^2))])+ ArcSin[(x/a)])],
x>a,(C/(2(((-(a^2)+(x^2)))^(1/4))))Exp[-(n+(1/2))((x/a) (Sqrt[((x^2)/(a^2))-1])- ArcCosh[(x/a)])]]]

Here is a graph showing the approximation of the wave function for n=0 and n=10. We see the divergence at the turning points.

pl1=Plot[{ψThHO[0,x],ψHO[0,x]},{x,-3,3},PlotStyle-> {Automatic,Dashing[0.02]},PlotRange-> {Automatic,{0,1.2}}];
pl2=Plot[{ψThHO[10,x],ψHO[10,x]},{x,-6,6},PlotStyle-> {Automatic,Dashing[0.02]},PlotRange-> {Automatic,{-0.5,0.8}}];
GraphicsRow[{pl1,pl2},ImageSize-> 7 72]

Our first task is to write a function which can compute levels and wave functions in WKB approximation numerically. We limit ourselves to simple well-type potentials, in which the search of turning points is simple, and for simplicity we assume that one turning point is for x< 0 while the second is positive. In other cases the reader can easily modify the function getTurningPoints below:

getTurningPoints[ene_,V_]:= Module[{aL,aR},
aR=x/.FindRoot[ene - V[x],{x,10}];
aL=x/.FindRoot[ene - V[x],{x,-10}];
{aL,aR}];

Now we need a function to compute action integrals

calcJ[ene_,V_]:= Module[{aR,aL,Jval},
{aL,aR}= getTurningPoints[ene,V];
Jval =(1/π)Re[NIntegrate[(Sqrt[2(ene - V[x])]),{x,aL,aR}]]];

and finally a function which implement quantization condition:

quantize[ene_?NumberQ,V_,n_]:=calcJ[ene,V] - (1/2) -n;

The WKB energy value will be given by something like:

getEneWKB[n_,V_]:=  y/. FindRoot[quantize[y,V,n] ==0,{y,0,0,200}]

The limits 0, 200 can be freely changed.

Let us check the procedure in the example of harmonic oscillator :

getEneWKB[#,V]&/@Range[0,5]

For the wave function we have just to write its definition. We suppose that energy has been computed. To keep functions independent we compute again turning points in this function.

eneWKB= getEneWKB[0,V]

calcPhase[ene_,V_]:=  Module[{aL,aR,T,C,σ,u},
{aL,aR} = getTurningPoints[ene,V];
σ =u/. First[NDSolve[{u'[z] ==  Sign[aR-z](Sqrt[2Abs[(-V[z]+ene )]]),u[aL]== 0},u,{z,5aL, 5aR}]];
T =  Re[ NIntegrate[(2/(Sqrt[2(ene - V[z])])),{z,aL,aR}]];
{{aL,aR,T},σ}]

ψWKB[x_,ene_,{aL_,aR_,T_},σ_]:= Module[{C},
C = 2/(Sqrt[T]);
Which[ x< aL,
(C/(2(((2(V[x]-ene)))^(1/4))))Exp[ σ[x]],
aL < x < aR,
(C/(((2(-V[x]+ene)))^(1/4)))Cos[σ[x] - (π/4)],
aR < x,
(C/(2(((2(V[x]-ene)))^(1/4))))Exp[ σ[x]- σ[aR]]]]

{{aL,aR,T},σ}=calcPhase[eneWKB,V];

Plot[{ψWKB[x,eneWKB,{aL,aR,T},σ],ψHO[0,x]},{x,-3,3.1},PlotStyle-> {Automatic,Dashing[0.02]},PlotRange-> {All,{0,1.2}}]

Anharmonic homogeneous oscillator

(x^4) case

As a second simple example we take

V[x_]:= (1/2)g (x^4)

In the problems solved for this chapter it is shown that WKB quantization condition leads to

eneHom[n_,g_]:= (1/2)(g^(1/3))(((n+(1/2)))^(4/3))(π^(2/3)) ((((( Gamma[(7/4)])/(Gamma[(5/4)]))))^(4/3))

First let us check our numerical procedure with these exact results

g=1.2;

getEneWKB[#,V]&/@Range[0,5]

eneHom[Range[0,5],g]

The results coincide.

Now we make a comparison with numerical computed eigenvalues. We use a variational procedure, fully explained in the chapter in variational calculations.

We compute matrices for operators (x^2), (x^4)etc. in a vector space of dimension dim = 200.

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

Then we solve the eigenvalue problem

{eigH,eigV} = Eigensystem[(1/2)p2 + (1/2)g x4];
order = Ordering[eigH];
eigH = eigH[[order]];
eigV= eigV[[order]];

eigH[[Range[6]]]

The agreement is satisfactory. We plot δE/E as a function of n

index = Range[0,50];

We remember that vectors in Mathematica use indexes starting from 1:

ListPlot[Transpose[{index, (eigH[[index+1]]-eneHom[index,g])/eigH[[index+1]]}]]

The classical periods of motion are easily computed, see problems

λ:= (g^(-1/6));
 ϵ0 := 1/2 /(λ^2)

classicalT[n_,g_]:= (2/(ϵ0))(((((eneHom[n,g])/(ϵ0))))^(-1/4))(((Sqrt[π]) Gamma[1+(1/4)])/(Gamma[(1/2)+(1/4)]))

and we can check the validity of the semiclassical relation

((∂(Subscript[E, n]))/(∂n))= ((2π)/T) ℏ

EofN = Interpolation[Transpose[{index,eigH[[index+1]]}]];

Plot[EofN'[x],{x,0,Max[index]},PlotRange-> All,AxesOrigin-> {0,1},
AxesLabel-> {"n","∂ E/∂ n"}, Epilog-> Point/@Transpose[{index,2π/classicalT[index,g]}]]

Wave functions

Let us compare wave functions. The numerical wave function can be obtained using the function getPsi on a grid:

state = 2;

xgrid = Range[-3.,3.,0.02];
Off[CompiledFunction::"cfn"]
psi = getPsi[eigV[[state+1]],xgrid];
On[CompiledFunction::"cfn"]

The WKB approximation is obtained as in the previous case

eneWKB=getEneWKB[state,V];

{{aL,aR,T},σ}=calcPhase[eneWKB,V];

listψWKB = ψWKB[#,eneWKB,{aL,aR,T},σ]&/@xgrid;

ListPlot[{Transpose[{xgrid,psi}],Transpose[{xgrid,listψWKB}]},Joined-> True,PlotStyle-> {Automatic,Dashing[0.02]},
PlotRange-> {All,{1.2Min[psi],1.2Max[psi]}}]

We see an excellent agreement except near the classical turning points, where we know that the approximation should fail.

(x^N)case

In the problems it is shown that for a generic potential 1/2 g (x^N)the Schrödinger equation can be reduced to the form

-((d^2)/(dz^2))ψ + (z^N)ψ = ϵ ψ

and energy eigenvalues can be written as

(Subscript[E, n]) = ((ℏ^2)/(2 m (λ^2))) (Subscript[ϵ, n]) ;  λ = (((((ℏ^2)/(m g))))^(1/(N+2)))

So what really matters are ϵ eigenvalues. In WKB approximation

(Subscript[ϵ, n]) = ((((n+(1/2)) (Sqrt[π])((Γ[(3/2)+(1/N)])/(Γ[1+(1/N)]))))^((2N)/(2+N))).

Let us call deg the degree of V[x], we can write

fϵ[n_]:= N@((((n+(1/2)) (Sqrt[π])((Gamma[(3/2)+(1/deg)])/(Gamma[1+(1/deg)]))))^((2deg)/(2+deg)))

For deg = 4 we have previous case. Let us compute again the eigenvalues. We compute matrices for dimension higher than dim and use multiplication of matrices to get (x^6)operator which will be used in a check below

deg = 4;

dim = 400;
{p2,x2,x4,x1}=N/@hermiteMatrices[dim+6];
x6 = x4.x2;
{p2,x2,x4,x1,x6}=Take[#,1;;dim,1;;dim]&/@{p2,x2,x4,x1,x6};

eigH = Sort[Eigenvalues[p2 + x4]];

index = Range[0,50];
index5 = Range[0,5];

δEoverE = (eigH[[index+1]]- fϵ[index])/eigH[[index+1]];

This is the table shown in the text

TableForm[Transpose[{eigH[[index5+1]],fϵ[index5],δEoverE[[index5+1]]}],
TableHeadings-> {index5,{"(Subscript[ϵ, n])","(Subscript[ϵ, n]^WKB)","δϵ/ϵ"}},TableSpacing-> {3,3}]

The same check can be done for others degree, as an instance deg = 6.

deg = 6;

eigH = Sort[Eigenvalues[p2 + x6]];

δEoverE = (eigH[[index+1]]- fϵ[index])/eigH[[index+1]];

TableForm[Transpose[{eigH[[index5+1]],fϵ[index5],δEoverE[[index5+1]]}],
TableHeadings-> {index5,{"(Subscript[ϵ, n])","(Subscript[ϵ, n]^WKB)","δϵ/ϵ"}},TableSpacing-> {3,3}]

Anharmonic oscillator

Analytic and numeric checks

In the problems of this chapter it has been shown that the most general parity invariant potential for an anharmonic oscillator with positive quadratic term is

V[x_]:= (1/2)(x^2) + (1/2)g (x^4)

It has also been shown that adiabatic invariant J and period can be computed in terms of elliptic integrals

J[E] = ((2 (Sqrt[g]))/(π ℏ))((Sqrt[(a^2)+(b^2)])/3)((b^2)EllipticK[((a^2)/((a^2)+(b^2)))] - ((b^2)-(a^2))EllipticE[((a^2)/((a^2)+(b^2)))]);

T[E] = (4/(Sqrt[g]))(1/(Sqrt[(a^2)+(b^2)]))EllipticK[((a^2)/((a^2)+(b^2)))].

EllipticE and EllipticK stand for elliptic integrals E and K respectively. a and b are given by

(a^2)=(( -1 + (Sqrt[1+ 8 g E]))/(2g)) ;  (b^2) = (( 1 + (Sqrt[1+ 8 g E]))/(2g))

and ±a are the classical inversion points. The classical momentum can be written as

p[x] = (Sqrt[2 E - (x^2)- (gx^4)]) = (Sqrt[g])(Sqrt[((a^2)-(x^2))((b^2)+(x^2))])

Our first task is to check these analytical expressions with our numerical routines. For the analytical expressions we have:

fJx4[ene_,g_]:= Module[{a,b},
a = (Sqrt[(( -1 + (Sqrt[1+ 8 g ene]))/(2g))]);b = (Sqrt[(( 1 + (Sqrt[1+ 8 g ene]))/(2g))]);
((2 (Sqrt[g]))/(π ))((Sqrt[(a^2)+(b^2)])/3)((b^2)EllipticK[((a^2)/((a^2)+(b^2)))] - ((b^2)-(a^2))EllipticE[((a^2)/((a^2)+(b^2)))])]

fTx4[ene_,g_]:= Module[{a,b},
a = (Sqrt[(( -1 + (Sqrt[1+ 8 g ene]))/(2g))]);b = (Sqrt[(( 1 + (Sqrt[1+ 8 g ene]))/(2g))]);
(4/(Sqrt[g]))(1/(Sqrt[(a^2)+(b^2)]))EllipticK[((a^2)/((a^2)+(b^2)))]]

Then we can write a routine for computing energy levels without performing numerical integration. The coupling constant now enters as a variable, instead of the whole potential

quantizex4[ene_?NumberQ,g_,n_]:=fJx4[ene,g]- (1/2) -n;
getEneWKBx4[n_,g_]:=  y/. FindRoot[quantizex4[y,g,n] ==0,{y,0,0,200},AccuracyGoal-> 6]

Let us compare the "analytical version" with the purely numerical one:

g = 0.7;

enex4=getEneWKBx4[#1,g]&/@Range[0,5]

getEneWKB[#1,V]&/@Range[0,5]

There is agreement. For the period the new version is

fTx4[#1,g]&/@enex4

while the numeric version is easily written:

calcT[ene_,V_]:= Module[{aR,aL,Tval},
{aL,aR}= getTurningPoints[ene,V];
Tval =2Re[NIntegrate[1/(Sqrt[2(ene - V[x])]),{x,aL,aR}]]];

calcT[#1,V]&/@enex4

WKB and variational results

We can perform the same checks done for the homogeneous potential. First let us compute eigenvalues and eigenvectors in the variational scheme. For s

g=0.2;

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

Then we solve the eigenvalue problem

{eigH,eigV} = Eigensystem[(1/2)p2 +(1/2)x2+ (1/2)g x4];
order = Ordering[eigH];
eigH = eigH[[order]];
eigV= eigV[[order]];

Then we compute WKB approximation up to level n = 50

index = Range[50];
enex4=getEneWKBx4[#1,g]&/@index;

and plot the relative error of the results :

ListPlot[Transpose[{index, (eigH[[index+1]]-enex4)/eigH[[index+1]]}],
AxesLabel-> {"n","δE/E"}]

We can also verify correspondence principle

EofN = Interpolation[Transpose[{index,eigH[[index+1]]}]];

classicalPeriods= fTx4[#1,g]&/@enex4;

Plot[EofN'[x],{x,0,Max[index]},PlotRange-> All,AxesOrigin-> {0,1},
AxesLabel-> {"n","∂ E/∂ n"}, Epilog-> Point/@Transpose[{index,2π/classicalPeriods}]]

Wave functions

As a last check we compare wave functions. We use function getPsi to compute the variational state

state = 3;

xgrid = Range[-3.,3.,0.02];
Off[CompiledFunction::"cfn"]
psi = getPsi[eigV[[state+1]],xgrid];
On[CompiledFunction::"cfn"]

The WKB approximation is obtained as in the previous case

eneWKB=getEneWKBx4[state,g];

{{aL,aR,T},σ}=calcPhase[eneWKB,V];

listψWKB = ψWKB[#,eneWKB,{aL,aR,T},σ]&/@xgrid;

We match the signs of the wave functions :

sign = Sign[listψWKB[[1]]]*Sign[psi[[1]]];

ListPlot[{Transpose[{xgrid,psi}],Transpose[{xgrid,sign*listψWKB}]},Joined-> True,PlotStyle-> {Automatic,Dashing[0.02]},
PlotRange-> {All,{1.2Min[psi],1.2Max[psi]}}]

There is a good agreement except near the turning points, again.