Airy functions and asymptotic (=WKB) forms

Aip[x_]:=1/(2 Sqrt[Pi Sqrt[x]]) Exp[-(2/3) x^1.5]

Aip2[x_] := 1/(Sqrt[Pi Sqrt[-x]]) Cos[(2/3)(-x)^1.5 - Pi/4]

airyplot=Plot[{AiryAi[x]},{x,-5,5},PlotRange->{{-5,5},{-.5,0.6}},PlotStyle->{{Thickness[0.01]},{Thickness[0.01],Dashing[0.1]}},AxesLabel->{"z"},PlotLabel->"Airy fcn vs. asymp. form"];

airywkb1plot=Plot[{Aip[x]},{x,0.1,5},PlotRange->{{-5,5},{-.5,0.6}},PlotStyle->{{Thickness[0.01],Dashing[0.02],Red}},AxesLabel->{"z"},PlotLabel->"Airy fcn vs. asymp. form"];

airywkb2plot=Plot[{Aip2[x]},{x,-5,-0.1},PlotRange->{{-5,5},{-.5,0.6}},PlotStyle->{{Thickness[0.01],Dashing[0.02],Red}},AxesLabel->{"z"},PlotLabel->"Airy fcn vs. asymp. form"];

Airy Plots

Show[airyplot,airywkb1plot,airywkb2plot]

Plot[{AiryAi[x],AiryBi[x]},{x,-10,5},PlotRange->{{-10,5},{-0.5,1}},PlotStyle->{{Thickness[0.01]}},PlotLabel->"Ai(z) vs. Bi(z)"]

SHO with hard wall

WKB result

Potential is (m omega^2/2)x^2 for |x|<a and infinite otherwise. Energies are measured in units of hbar omega/2, lengths in units of sqrt(hbar/(m omega)), which is Sqrt[2] times the spread of the g.s. wavefunction in a SHO.

theta0[e_,a_]= ArcSin[a/Sqrt[e]];
lhs[e_,a_]=(theta0[e,a]+Sin[2theta0[e,a]]/2) (e/Pi);

Module to find n'th WKB energy for width a:

rootwkb[n_,a_]:= Module[{x,y}, y=FindRoot[lhs[x,a]==n,{x,5.1}];x/.y]

Results for a=1. Compare to SHO results which are just 1, 3, 5, 7

Table[{n,rootwkb[n,1]},{n,1,5}]

If widen the well, one eventually can find no oscillating solution for n=1. (To obtain this in WKB one would have to match to a decaying part outside.) The potential at the wall is given by a^2. Thus, for a=2, all energies
at or below E=4 have classical turning points inside the walls, and cannot be obtained without matching to decaying exponentials. Thus the following energies are only valide for n=3,4 & 5. Note that the energies of the these three states are reduced substantially compared to a=1. This is because these states are much less "squeezed" by the walls.

Table[{n,rootwkb[n,2]},{n,1,5}]

Plotting WKB wavefunction:

wkbwf[x_,e_,a_]:=Module[{thetaf,thetai,phase},thetaf =ArcSin[x/Sqrt[e]];thetai=theta0[e,a];phase=(thetaf+thetai)/2+Sin[2 thetaf]/4 +Sin[2 thetai]/4;Sin[e phase]/Sqrt[e-a^2]]

plotwkb[n_,a_]:=Module[{e},e=rootwkb[n,a];Plot[wkbwf[x,e,a],{x,-a,a},PlotStyle->{Red,Thickness[0.01]}]]

plotwkb[4,1]

V=0 results

Energies if V=0

freee[n_,a_]:=N[ (n/a)^2 (Pi/2)^2]

Table[freee[n,1],{n,1,5}]

Exact solutions

Exact solutions (using trick of treating energy as variable to be solved for). Input is starting guess for energy. Output is plot of wavefunction with label giving exact energy. Use WKB values as initial guesses.

exact[e0_,a_]:=Module[{u,x,e,ee},sol=NDSolve[{u''[x]==u[x](x^2-e[x]),e'[x]==0,u'[-a]==1,u[-a]==0,u[a]==0},{u,e},{x,-a,a},Method->{"Shooting","StartingInitialConditions"->{e[0]==e0,u[0]==1,u'[0]==1}}];ee=e[0]/.First[sol];Plot[u[x]/.First[sol],{x,-a,a},PlotStyle->{{Thickness[0.01]}},PlotLabel->Evaluate[ee]]]

eexact[n_,a_]:=Module[{u,x,e,ee},sol=NDSolve[{u''[x]==u[x](x^2-e[x]),e'[x]==0,u'[-a]==1,u[-a]==0,u[a]==0},{u,e},{x,-a,a},Method->{"Shooting","StartingInitialConditions"->{e[0]==rootwkb[n,a],u[0]==1,u'[0]==1}}];ee=e[0]/.First[sol]]

WKB vs. exact energies for first 4 states with a=1

Table[{rootwkb[n,1],eexact[n,1]},{n,1,4}]

Plots

What follows are exact wavefunctions for first 4 excited states(blue) along with WKB wavefunctions (red)

exact1=exact[2,1];Show[plotwkb[1,1],exact1]

exact2=exact[10,1];Show[exact2,plotwkb[2,1]]

exact3=exact[22,1];Show[exact3,plotwkb[3,1]]

exact4=exact[40,1];Show[exact4,plotwkb[4,1]]

Show[exact2,exact1,exact3,exact4,PlotRange->{{-1,1},{-.4,.7}}]

Here's what the exact g.s. looks like for the walls twice as far apart. Note that the energy is almost down to 1 (i.e. h-bar omega/2), although the wavefunction is clearly quite distorted from the oscillator ground state.

exact[2,2]

The distortion of the first excited state is more substantial (with the energy being
increased from 3 -> 3.53).

exact[4,2]

Here is the case of a much wider potential, where the wavefunction and energy have basically collapsed to the SHO values.

exact[1,4]

Quartic potential

Potential is V_0 (x/x_0)^4. Energies measured in units of V_0 (L/x_0)^4, where L^3=h-bar x_0^2/sqrt(2mV_0).

WKB

phaseq =NIntegrate[Sqrt[1-y^4],{y,-1,1}]

rootwkbq[n_]:=N[ ( (n+1/2) Pi/phaseq)^(4/3)]

WKB energies:

Table[rootwkbq[n],{n,0,5}]

uwkbq[x_,e_]:=Module[{y0},y0=x/e^(1/4);(1/(1-y0^4)^(1/4))Cos[e^(3/4)NIntegrate[Sqrt[1-y^4],{y,-1,y0}]-Pi/4]]

Plotting WKB wavefunctions, normalized to unity at position x0:

plotq[n_,x0_]:=Module[{e,eps,u0},eps=0.01;e=rootwkbq[n];u0=uwkbq[x0,e];Plot[uwkbq[y,e]/u0,{y,-e^(1/4)+eps,e^(1/4)-eps},PlotStyle->{{Red,Thickness[0.01]}}]]

Can see the impact of the turning points:

plotq0=plotq[0,0];

plotq1=plotq[1,0.05];

plotq2=plotq[2,0];

plotq3=plotq[3,0.05];

plotq4=plotq[4,0];

plotq8=plotq[8,0];

WKB wavefunctions in E<V region:

uwkbqout[x_,e_]:=Module[{y0},y0=x/e^(1/4);(1/2)(1/(y0^4-1)^(1/4))Exp[-e^(3/4)NIntegrate[Sqrt[y^4-1],{y,1,y0}]]]

plotqout[n_,rat_,x0_]:=Module[{e,eps,u0},eps=0.002;e=rootwkbq[n];u0=uwkbq[x0,e];Plot[uwkbqout[y,e]/u0,{y,e^(1/4)+eps,rat e^(1/4)},PlotStyle->{{Red,Thickness[0.01]}}]]

plotqout0=plotqout[0,3,0];

plotqout1=plotqout[1,2.4,-0.05];

plotqout2=plotqout[2,2,0];

plotqout3=plotqout[3,1.9,-0.05];

plotqout4=plotqout[4,1.8,0];

plotqout8 = plotqout[8,1.5,0];

Exact energies and wavefunctions

Solving exactly. "rat" is a fudge factor determining how far out one should set the wavefunction to zero compared to the turning point. Could be more sophisticated using the WKB wavefunction itself.

First for even parity eigenvalues:

eexactevenq[n_,rat_]:=Module[{u,x,e,ee,ep,xp,sol},ep=rootwkbq[n];xp=rat ep^(1/4);sol=NDSolve[{u''[x]==u[x](x^4-e[x]),e'[x]==0,u'[0]==0,u[-xp]==0,u[xp]==0},{u,e},{x,-xp,xp},Method->{"Shooting","StartingInitialConditions"->{e[0]==ep,u[0]==1,u'[0]==0}}];ee=e[0]/.First[sol]]

Here are the lowest three (note that n=0 is allowed):

eexact0=eexactevenq[0,3];eexact2=eexactevenq[2,2]; eexact4=eexactevenq[4,1.8];

{eexactevenq[0,3],eexactevenq[2,2],eexactevenq[4,1.8]}

Compare these to the wkb values:

Table[rootwkbq[n],{n,0,4,2}]

Next the parity odd states (enforced by u[0]=0; note also the initial conditions are changed):

eexactoddq[n_,rat_]:=Module[{u,x,e,ee,ep,xp,sol},ep=rootwkbq[n];xp=rat ep^(1/4);sol=NDSolve[{u''[x]==u[x](x^4-e[x]),e'[x]==0,u[0]==0,u[-xp]==0,u[xp]==0},{u,e},{x,-xp,xp},Method->{"Shooting","StartingInitialConditions"->{e[0]==ep,u[0]==0,u'[0]==1}}];ee=e[0]/.First[sol]]

eexact1=eexactoddq[1,2.4];eexact3=eexactoddq[3,1.9];

{eexactoddq[1,2.4],eexactoddq[3,1.9]}

{rootwkbq[1],rootwkbq[3]}

Plotting routines for exact wavefunctions

plotexactevenq[n_,rat_]:=Module[{u,x,e,ee,ep,xp,sol,u0},ep=rootwkbq[n];xp=rat ep^(1/4);sol=NDSolve[{u''[x]==u[x](x^4-e[x]),e'[x]==0,u'[0]==0,u[-xp]==0,u[xp]==0},{u,e},{x,-xp,xp},Method->{"Shooting","StartingInitialConditions"->{e[0]==ep,u[0]==1,u'[0]==0}}];ee=e[0]/.First[sol];u0=u[0]/.First[sol];Plot[u[x]/u0/.First[sol],{x,-xp,xp},PlotStyle->{{Thickness[0.01]}}]]

plotexactoddq[n_,rat_]:=Module[{u,x,e,ee,ep,xp,sol,u0},ep=rootwkbq[n];xp=rat ep^(1/4);sol=NDSolve[{u''[x]==u[x](x^4-e[x]),e'[x]==0,u[0]==0,u[-xp]==0,u[xp]==0},{u,e},{x,-xp,xp},Method->{"Shooting","StartingInitialConditions"->{e[0]==ep,u[0]==0,u'[0]==1}}];ee=e[0]/.First[sol];u0=u[.05]/.First[sol];Plot[u[x]/u0/.First[sol],{x,-xp,xp},PlotStyle->{{Thickness[0.01]}}]]

plotexactq0=plotexactevenq[0,3];

plotexactq1=plotexactoddq[1,2.4];

plotexactq2=plotexactevenq[2,2];

plotexactq3=plotexactoddq[3,1.9];

plotexactq4=plotexactevenq[4,1.8];
plotexactq8 = plotexactevenq[8,1.5];

Plotting Airy function in connection region

airy[x_,e_,normpt_]:=Module[{xt,alpha},xt=e^(1/4);alpha=(0.25)^(1/3);AiryAi[(x-xt)/alpha]/AiryAi[(normpt-xt)/alpha]]

plotairyeven[n_,rat_]:= Module[{xt,alpha,xp,u,u0,ut,x,e,ee,ep,sol},ep=rootwkbq[n];xp=rat ep^(1/4);sol=NDSolve[{u''[x]==u[x](x^4-e[x]),e'[x]==0,u'[0]==0,u[-xp]==0,u[xp]==0},{u,e},{x,-xp,xp},Method->{"Shooting","StartingInitialConditions"->{e[0]==ep,u[0]==1,u'[0]==0}}];e=e[0]/.First[sol];xt=e^(1/4);u0=u[0]/.First[sol];ut=u[xt]/.First[sol];xp=rat xt;alpha=(0.25)^(1/3)/xt;ai0=AiryAi[0];Plot[AiryAi[(x-xt)/alpha]ut/(u0 ai0),{x,0,xp},PlotStyle->{{Orange,Thickness[0.01]}}]]

plotai0=plotairyeven[0,3];

plotai2=plotairyeven[2,2];

plotai4=plotairyeven[4,1.8];

plotai8=plotairyeven[8,1.5];

Plots comparing exact, Airy and WKB wavefunctions

The following plots show the "exact" wavefunction in blue, the inside and outside WKB solutions in red and the Airy function in orange (the latter only in some plots). The WKB curves are normalized to match the normalization of the exact solution either at 0 (even cases) or 0.05 (odd cases). The relative normalization of the inside and outside  WKB solutions are given by the connection formula. The Airy function is normalized to match the exact wavefunction at the turning point, but is otherwise fixed.

One sees that the matching is very poor for n=0, but gets progressively better as n increases. Nevertheless, even for n=4 there is not a region where the inner WKB solution matches onto the Airy function. This appears to be mostly a problem with the overall amplitude, since the WKB solution matches the phase of the exact solution well.

Show[plotexactq0,plotq0,plotqout0,plotai0]

Show[plotexactq1,plotq1,plotqout1]

Show[plotexactq2,plotq2,plotqout2,plotai2]

Show[plotexactq3,plotq3,plotqout3]

Show[plotexactq4,plotq4,plotqout4,plotai4]

Show[plotexactq8,plotq8,plotqout8,plotai8]