WKB approximation

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study the double well anharmonic oscillator in the WKB approximation, compare the results with variational computations. Study in particular the small coupling limit.

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

Introduction

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}]];

getscale[n_,g_,μ_]:= x/.FindRoot[(2n+1)(x-(μ/(x^3)))-((3g)/(x^5))(2(n^2)+2n+1),{x,(((3g))^(1/3))}]

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 model and the WKB approximation

This is the potential to be studied

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

Block[ {μ = -1,  g=0.1,ene=-0.5,a,b},
a = (Sqrt[((1-(Sqrt[1+8 ene g]))/(2g))]);b = (Sqrt[((1+(Sqrt[1+8 ene g]))/(2g))]);
 Plot[V[x],{x,-1.5(Sqrt[((-μ)/g)]),1.5(Sqrt[((-μ)/g)])},AxesLabel-> {"x","V[x]"},PlotRange-> {{-5,5},{-1.5,2}},Epilog-> {Dashing[0.015],Red,Line[{{-b,ene},{-a,ene}}],Line[{{b,ene},{a,ene}}],
Black,Style[{Text["b",{b+0.4,ene}],Text["a",{a-0.4,ene}]},Medium]} ]]

The transition points a and b are given by

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

In problem [9] an analytic expression for action integral, classical period and penetration factor was given

J = ((Sqrt[g])/π)(∫_a^b)(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = ((Sqrt[g])/π)(b/3)(((a^2)+(b^2))EllipticE[(((b^2)-(a^2))/(b^2))] - 2 (a^2)EllipticK[(((b^2)-(a^2))/(b^2))]);

T = (2/(Sqrt[g]))(∫_a^b)(1/(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]))dx = (2/(Sqrt[g]))(1/b)EllipticK[(((b^2)-(a^2))/(b^2))];

K = Exp[-2 D];  D = (Sqrt[g])(∫_0^a)(Sqrt[( (a^2)- x^2)((b^2)-(x^2))])dx = (Sqrt[g]) (b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] -  ((b^2) - (a^2))EllipticK[((a^2)/(b^2))])

J fixes the quantization condition inside the well :

J = n+ 1/2

while K and T give the energy splitting between ground state and first excited state

ΔE =  ((2 ℏ)/T)  Exp[-(1/ℏ)(∫_(-a)^(+a))|p[x]|dx]≡ ((2 ℏ)/T) K .

Here is the Mathematica translation of the above formulas:

fJx4[ene_,g_]:= Module[{a,b},
a = (Sqrt[(( 1 - (Sqrt[1+ 8 g ene]))/(2g))]);b = (Sqrt[(( 1 + (Sqrt[1+ 8 g ene]))/(2g))]);
Re[(( (Sqrt[g]))/(π ))(b/3)(((b^2)+(a^2))EllipticE[(((b^2)-(a^2))/(b^2))] - 2(a^2)EllipticK[(((b^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))]);
 (2/(Sqrt[g]))(1/b)EllipticK[(((b^2)-(a^2))/(b^2))]];

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

Checks

Our first task is to provide some numerical check for the above expressions. We can follow the procedure used in the study of usual anharmonic oscillator. We just report here functions used in notebook [], simply redefine the function getTurningPoints by using their analytical expression. We leave intact the syntax but use g as a variable, instead of V[x]. Here aL and aR are turning points inside a well:

getTurningPoints[ene_,g_]:= Module[{aL,aR},
aR=(Sqrt[(( 1 + (Sqrt[1+ 8 g ene]))/(2g))]);aL= (Sqrt[(( 1 - (Sqrt[1+ 8 g ene]))/(2g))]);{aL,aR}];

Here we make a direct computation of integrals:

calcJ[ene_,g_]:= Module[{aR,aL,Jval},
{aL,aR}= getTurningPoints[ene,g];
Jval =(1/π)Re[NIntegrate[(Sqrt[2(ene -(((-1)/2)(x^2)+ (g/2)(x^4)))]),{x,aL,aR}]]];

calcK[ene_,g_]:= Module[{aL,Jval,dF},
aL= (Sqrt[(( 1 - (Sqrt[1+ 8 g ene]))/(2g))]);
dF =Re[NIntegrate[(Sqrt[2(-ene +(((-1)/2)(x^2)+ (g/2)(x^4)))]),{x,0,aL}]];
Exp[-2 dF]];

calcT[ene_,g_]:= Module[{aR,aL,Tval},
{aL,aR}= getTurningPoints[ene,g];
Tval =2Re[NIntegrate[1/(Sqrt[2(ene -(((-1)/2)(x^2)+ (g/2)(x^4)))]),{x,aL,aR}]]];

We prepare a bunch of random couples (E,g) to perform the check. E is chosen between 0 and the minimum of the potential, - 1/8g.

SeedRandom[7];
gRan = RandomReal[{0.01,0.2},{30}];
eneRan= - RandomReal[{0.01,1},{30}]/(8 gRan);
list=Transpose[{eneRan,gRan}];

Max[(fJx4@@ #1&/@list) -( calcJ@@ #1&/@list)]

Max[(fKx4@@ #1&/@list) -( calcK@@ #1&/@list)]

Max[(fTx4@@ #1&/@list) -( calcT@@ #1&/@list)]

The checks are satisfactory.

WKB and variational results

listg = Range[0.02,0.18,0.01];

dim=200;

{p2,x2,x4} = N/@hermiteMatrices[dim,0];
eneEven = Table[g1= listg[[k]];α = getscale[dim,g1,-1];
Sort[Eigenvalues[((α^2)/2)p2 - (1/(2(α^2)))x2+ (1/(2(α^4)))g1 x4,-50,Method-> "Arnoldi"]][[1]],{k,1,Length[listg]}];

{p2,x2,x4} = N/@hermiteMatrices[dim,1];
eneOdd = Table[g1= listg[[k]];α = getscale[dim,g1,-1];
Sort[Eigenvalues[((α^2)/2)p2 - (1/(2(α^2)))x2+ (1/(2(α^4)))g1 x4,-50,Method-> "Arnoldi"]][[1]],{k,1,Length[listg]}];

δE = -eneEven+eneOdd;

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

eneWKB=getEneWKBx4[0,#1]&/@listg;

listEneg = Transpose[{eneWKB,listg}];
classicalPeriods = fTx4@@#1&/@listEneg;
δEWKB = 2/classicalPeriods * (fKx4@@#1&/@listEneg);

ListPlot[Transpose[{listg,δEWKB}],Joined-> True,
AxesLabel-> {"g","δE"},
Epilog-> Point/@Transpose[{listg,δE}]]

TableForm[Transpose[{δE,δEWKB}],TableHeadings-> {listg,{"δE","δEWKB"}},
TableSpacing-> {1.5,3}]

The small coupling limit

The analytical results used in this section are derived in problem [11] to which the interested reader is referred.

The analytical results

For small g the barrier becomes very high and the penetration coefficient K very small. This regime can be analyzed by perturbation theory and WKB. The main results are the following.

In this model the real expansion parameter, at fixed quantum number, is ℏ g, i.e. small ℏ leads to small g expansion.

In usual WKB the prefactor in the formula ΔE ∝ K is essentially unknown. Using uniform WKB approximation it can be computed, in the limit →0

(Subscript[ΔE, n]) = 2 ((((π/2)))^(1/2))((((((n+1/2))/e)))^(n+(1/2)))(1/(n!)) ((ℏω)/π)K

In perturbation theory (Subscript[ΔE, 0]) for the ground state can be computed at arbitrary order in perturbation theory, the first orders in the expansion being

(Subscript[ΔE, 0]) =(1/(Sqrt[π]))(2^(9/4))((e^(-((Sqrt[2])/(3g))))/(Sqrt[g]))(1-((71 g)/(24 (Sqrt[2])))-((6299 (g^2))/2304))

The results () and () agree in lowest order perturbation theory.

In the following we will compare high precision numerical results for small g with predictions () and ().

The prefactor

The prefactor in () is always close to 1, this explain why the naive WKB formula gives a good agreement with numerical data.

prefactor[n_]:=  2 ((((π/2)))^(1/2))((((((n+1/2))/E)))^(n+(1/2)))(1/(n!))

ListPlot[prefactor[Range[0,20]],PlotRange-> {0.9,1.1},AxesLabel-> {"n","p[n]"}]

The largest correction is for the ground state, and is (Sqrt[π/e])∼1.07, asymptotically the prefactor goes to 1 and naive WKB formula is recovered.  The formula () has corrections of order g.  For relatively high g its reliability is questionable. As a matter of fact it appears that for these g the naive formula is more accurate, this seems accidental. The perturbative formula for the ground state is always accurate: we remember that in this problem we have splitting only for g < 0.2, i.e. the coupling cannot be "large" (in units where ℏ=1).

Low precision results

The previous observations are substantiated by making a table of previous results

dEGSpert[g_]:= (1/(Sqrt[π]))(2^(9/4))((E^(-((Sqrt[2])/(3g))))/(Sqrt[g]))(1-((71 g)/(24 (Sqrt[2])))-((6299 (g^2))/2304))

δEWKBCorr = prefactor[0]*δEWKB;

TableForm[Transpose[{δE,δEWKB,δEWKBCorr,dEGSpert[listg]}],TableHeadings-> {listg,{"δE","δEWKB","δEWKBCorrected","δEpert"}},
TableSpacing-> {1.5,3}]

The critial check is in the low g region, where we see that the prefactor improve the agreement.

High precision results

To compute ΔE numerically at small g we need an high precision on individual levels. This problem has already be faced in the chapter on variational calculations where the two following routines have been written to perform an high precision calculation of energy levels.

First we compute matrix elements with a precribed precision

computeEneHighPrecision[state_,parity_,α_,dim_,prec_:MachinePrecision]:= 
Module[{p2,x2,x4,ene0},
{p2,x2,x4}=N[hermiteMatrices[dim,parity],prec];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
ene0=Sort[Eigenvalues[Normal[Ham]]][[state]]];

If our sought precision is precisionGoal, it is reasonable to work with higher and higher dimensions up the precision sought is obtained. The following function simply make a loop, with a while control which stop execution when the precision is obtained or when the dimension go beyond the limit dimMax. An optional boolean argument printResults print intermediate results if settled to True (the default)

energyEstimate[state_,parity_,α_,precGoal_,dimMax_,printResults_:True]:=
Module[{dim,prec,ene,deltaE,enenew,precEne,enesafe,risOld ,deltaEold,precEneold},
dim = 30;
shiftdim = 10;
prec = precGoal+10;
precEne = prec;
ene = computeEneHighPrecision[state,parity,α,dim,prec];
enesafe = ene;
deltaE = 1;
While[ 
dim ≤ dimMax -shiftdim&& precEne ≥ precGoal && deltaE ≥ 10^(-precGoal),
dim = dim+shiftdim;
enenew = computeEneHighPrecision[state,parity,α,dim,prec];
deltaEold =  Norm[ene - enesafe];
precEneold = Precision[ene];
deltaE = Norm[ene - enenew];
precEne = Precision[enenew];
enesafe = ene;
ene = enenew;
risOld = {enesafe,dim-shiftdim,precEneold,deltaEold};
If[printResults== True,
Print[TableForm[
{{"dim=  " <> ToString[dim],
NumberForm[enenew,Ceiling[precGoal]+1],
" prec.= "<>ToString[precEne],"δE:",ScientificForm[deltaE,4]
}}]];];];
If[dim ≤ dimMax && deltaE < deltaEold,{ene,dim,precEne,deltaE},risOld]
]

This fixes our parameters

state=Range[1,20];
precGoal=50;
maxDim = 150;

g=5/1000;

dim=200;
precParameters=100;

We compute first the energies for even states

μ=-1;parity=0;
α = getscale[dim/2,g,μ];
{g,μ,α} = SetPrecision[#,precParameters]&/@{g,μ,α};

({ene0,lastDim,precEne,deltaE}=energyEstimate[state,parity,α,precGoal,maxDim,False]);

Then for odd states

μ=-1;parity=1;
α = getscale[dim/2,g,μ];
{g,μ,α} = SetPrecision[#,precParameters]&/@{g,μ,α};

({ene1,lastDim,precEne,deltaE}=energyEstimate[state,parity,α,precGoal,maxDim,False]);

ScientificForm[δE =N[ ene1-ene0]]

These are the WKB energies

ind = state[[1;;20]]-1;
eneWKB=getEneWKBx4[#1,g]&/@ind

and their relative WKB splitting

listg= ConstantArray[g,Length[ind]];

listEneg = Transpose[{eneWKB,listg}];
classicalPeriods = fTx4@@#1&/@listEneg;
δEWKB = 2/classicalPeriods * (fKx4@@#1&/@listEneg);

If we introduce the prefactor of equation () we get

δEWKBCorr = δEWKB*prefactor[state];

Here is a comparison of the results :

TableForm[Transpose[{δE,δEWKBCorr,δEWKB}],TableHeadings-> {state,{"δE","δEWKBCorrected","δEWKB","δEpert"}},
TableSpacing-> {1.5,3}]

We see a general agreement of the order of 2. (10^(-3))in relative precision, this is reasonable as g = 5. (10^(-3))and errors are of order 𝒪(g). Higher order corrections in g clearly improve the result for the ground state:

N@dEGSpert[g]

Let us consider another peculiarity of this kind of system. A strict perturbative expansion of eq.() gives

dEpert[n_,g_]:= (((2^((9/4)+((5 n)/2))) )/((Sqrt[π]) n!))(1/(g^(n+(1/2))))(E^(-((Sqrt[2])/(3 g))))

While this expression is formally correct its blind application can give bad results for higher states. In effect as it is evident even in perturbation theory the effective coupling in these cases is (n g). As an instance the perturbative expansion of energy levels for an anharmonic oscillator is

(Subscript[E, pert]) = ℏ(n+(1/2))+ (3/8) g (1+2 n+2 (n^2)) (ℏ^2)-(1/32) (g^2) (21+59 n+51 (n^2)+34 (n^3)) (ℏ^3)

Then for high states E ∼ (n+1/2) (1 + 𝒪(n g) ) and it is n g that must be small. This parameter appear in the exponent for the transition amplitude K, and even a small error can give rise to bad results for energy splitting. This is evident in the following table, where previous results are compared with those of last formula:

TableForm[Transpose[{δE,δEWKBCorr,δEWKB,N@dEpert[state-1,g]}],TableHeadings-> {state,{"δE","δEWKBCorrected","δEWKB","δEpertNaive"}},
TableSpacing-> {1.5,3}]