Accuracy of the Wentzel-Kramers-Brillouin (WKB) Quantization Condition for Symmetric Power Law Potentials

Frank J. Kampas
Physicist at Large Consulting

Frank @ Physicistatlarge.com

Abstract

It is well known that the WKB quantization condition gives the correct energy levels for the quantum mechanical simple harmonic oscillator (quadratic potential) and that it gives a good approximation in general for states with a large quantum number.  For symmetric power law potentials V(x) ∝ | ((x )^p) | , it is shown that the WKB energy levels  ∝ (((n + (1/2)))^((2p)/(2+p))) where n is the quantum number. Furthermore, the value of the WKB prediction divided by the actual value for the ground state energy decreases with increasing p, with values of approximately  1.09, 1.0, 0.90 ,0.82, 0.75,and 0.25 for p = 1, 2, 3, 4, 5, ∞.  The accuracy of the inverse WKB quantization condition for some of these potentials is also discussed.

Introduction

The WKB approximation is derived by writing the solution of the time-independent Schrodinger equation as A Exp( i S(x)/ℏ ), expanding the resulting differential equation in powers of ℏ and only going to first order (Schiff, p. 277).  The solution to the approximate equation is broken into 3 parts, an oscillatory region between the classical turning points (where the total energy is 0) and exponentially decaying regions outside the turning points.

According to the WKB approximation, the quantization condition for a one-dimensional system is given by

2(∫_(Subscript[x, 1])^(Subscript[x, 2]))((( 2m (  E-V(x) ) ))^(1/2))dx = ( n + (1/2) ) h

where  m is the particle mass, E is the total energy, V(x) is the potential function, (Subscript[x, 1]) and (Subscript[x, 2]) are the turning points, h is Planck’s constant and n is the number of nodes in the WKB wave function.

WKB Energy Levels for a Symmetric Polynomial Potential

For a potential which is symmetric in x so that V(x) = V(-x), the quantization condition can be rewritten as

(∫_0^(Subscript[x, 1]))(((  E-V(x) ))^(1/2))dx = ( n + (1/2) )(( h)/(4(Sqrt[2m])))

If the potential is a monotonic function of x for x > 0,  the variable of integration can be changed from x to V giving

(∫_0^E)(((  E-V  ))^(1/2))(dx/dV)dV = ( n + (1/2) )(( h)/(4(Sqrt[2m])))

Integrating by parts to get an expression containing x(V) rather than x’(V) gives

(∫_0^E)(((  E-V  ))^(-(1/2)))x(V)dV =  ( n + (1/2) )(( h)/(2(Sqrt[2m])))

replacing h with 2 π ℏ gives

(∫_0^E)(((  E-V  ))^(-(1/2)))x(V)dV =  ( n + (1/2) )(( π ℏ)/(Sqrt[2m]))

Assume that the potential function can be written as V(x) = k ((Abs(x))^p).  Then the integral  (∫_0^E)(((  E-V  ))^(-(1/2)))x(V)dV is

int = Integrate[(((e-v))^(-1/2))((((v/k)))^(1/p)),{v,0,e},Assumptions -> {e > 0, p > 0}]

Therefore the energy levels go as n + (1/2) raised to ((2 p)/(2+p)) power.

energies[p_] = e /. Solve[PowerExpand[int] == (n+(1/2))((π ℏ)/(Sqrt[2m])),e][[1]]//FullSimplify

Quantization Condition for the Dimensionless Problem

For potentials for which there is no analytic solution, it is useful to write the Schrodinger equation in dimensionless form.  For a symmetric power law potential, k (x^p), the equation is

Seq =  - ((ℏ^2)/(2m))ψ''[x] + (k (x^p)-e)ψ[x] == 0;

It is not necessary to write the potential as k ((Abs(x))^p)) if the equation is solved from x = 0 to x = ∞, rather than from -∞ to ∞.  In this approach, the solution will be even under inversion about x = 0 if ψ(x) = 1 and ψ’(x) = 0 or odd if ψ(x) = 0 and ψ’(x) = 1.

The scale change consists of replacing x with ρ y and e with σ ϵ.

Seq1 = Seq /. {x -> ρ y, ψ''[x] -> (ρ^(-2))ψ''[y], ψ[x] -> ψ[y], e -> σ ϵ}

The desired result is obtained if ρ is equal to ( ((ℏ^2)/(2k m)) )^(1/(2+p)) and  σ is equal to ( k  (  (ℏ/(Sqrt[2m])) )^p )^(2/(2+p))

Seq2 = Seq1 /. { ρ-> (((((ℏ^2)/(2 k m))))^(1/(2+p))), σ -> (((k ((((ℏ/(Sqrt[2 m]))))^p)))^(2/(2+p)))}//PowerExpand //FullSimplify

Therefore, the resulting equation is ψ’’(y) + (ϵ - (y^p)) ψ(y) = 0

The same scale change can be applied to the WKB quantization condition

((4(Sqrt[2 m]))/(2 π ℏ))(∫_0^(Subscript[x, 1]))(((   E-V(x) ))^(1/2))dx = ( n + (1/2) )

The integrand is

wkbint = ((4 (Sqrt[2m]))/(2 π ℏ))(Sqrt[e - k (x^p)])dx /. {x -> ρ y,dx -> ρ dy, e -> σ ϵ}/. { ρ-> (((((ℏ^2)/(2 k m))))^(1/(2+p))), σ -> (((k ((((ℏ/(Sqrt[2 m]))))^p)))^(2/(2+p)))}//PowerExpand//FullSimplify

Moving the term outside the square root (except for dy) inside the square root gives

((dy  (Sqrt[((2^(-(p/(2+p)))) (k^(2/(2+p))) (m^(-(p/(2+p)))) (-(y^p)+ϵ) (ℏ^((2 p)/(2+p))))((((2^((3/2)-(1/(2+p)))) (k^(-(1/(2+p)))) (m^(p/(4+2 p))) (ℏ^(-(p/(2+p))))))^2)]))/π)//PowerExpand//FullSimplify

Therefore

(∫_0^(Subscript[y, 1]))  (((ϵ - (y^p)))^(1/2))dy =(π/2)( n + (1/2) )

Performing the integral gives

scaledint = Integrate[(((ϵ - (y^p)))^(1/2)),{y,0,(ϵ^(1/p))},Assumptions -> {ϵ > 0, p > 0}]

which enables calculation of the scaled energies as a function of n

scaledenergies[p_] = ϵ /. Solve[scaledint == (π/2)(n+ (1/2)),ϵ][[1]]//FullSimplify

Note that the coefficient of (2n + ((1))^((2p)/(2+p))) is approximately 1 and decreases gradually with p

Table[N @ scaledenergies[p], {p,1,5}]

In the limit as p -> ∞, the coefficient goes to about 0.6

Limit[scaledenergies[p], p -> ∞]

N[Simplify @ ((%)/(((2n+1))^2))]

Analytic Solutions

In this section, the results of the WKB quantization condition for symmetric linear and quadratic potentials are compared to literature results for the energy levels.

Quadratic Potential (SHO)

For the SHO, with V(x) = (1/2)m (ω^2) (x^2), the result agrees with the result from solving the Schrodinger equation, (Subscript[E, n]) = (n + 1/2) ℏ ω,  as is well known.

energies[2]/. {k -> (1/2)m (ω^2)}//PowerExpand

As expected, the same result is obtained from the scaled energies, using the conversion factor from ϵ to e

(((k ((((ℏ/(Sqrt[2 m]))))^p)))^(2/(2+p)))scaledenergies[2]/. {k -> (m/2)(ω^2)}/. p -> 2//PowerExpand

Linear Potential

For a  symmetric linear potential V(x) = k |x], the WKB quantization condition gives

energies[1]//PowerExpand

Numerically, this is

N[%]

which can be written as

1.77068(((n + 1/2))^(2/3)) (  ((ℏ k)/(Sqrt[2m])) )^(2/3)

The Schrodinger equation for the linear symmetric potential be solved analytically. The energy levels are given by (Subscript[E, i]) =  ( ((ℏ k)/(Sqrt[2m])) )^(2/3) * (Subscript[y, i]) where Ai(-(Subscript[y, i])) = 0 or (Subscript[A, i])'(-(Subscript[y, i])) = 0 where Ai is one of the solutions to the Airy equation.(Robinett, p. 210).

The energy scaling factor for p = 1 is

(((k ((((ℏ/(Sqrt[2 m]))))^p)))^(2/(2+p)))/. p -> 1

which is equal to  ( ((ℏ k)/(Sqrt[2m])) )^(2/3).  Therefore the values of y from Ai(-(Subscript[y, i])) = 0 or (Subscript[A, i])'(-(Subscript[y, i])) = 0 can be compared directly to the scaled results.

Mathematica has a built-in function for the zeroes of the Airy function (AiryAiZero).  The zeroes of the derivative can be calculated as follows (from the documentation for AiryAiPrime)

b[n_]:=x/.FindRoot[AiryAiPrime[x],{x,-(((((3π)/8)(Max[0,4n-1])))^(2/3)),-(((((3π)/8)(4n+3)))^(2/3))}]

zs are the negative of the first 200 zeroes of Airy and AiryPrime.

zs = Sort[-Join[b[#]& /@ Range[0,99],N @ AiryAiZero[Range[100]]]];

For large values of n, it is known that (Subscript[E, n]) ∝ (n^(2/3))  (Robinett, p. 211), in agreement with an exponent of  ((2 p)/(2+p)) for p = 1.  Therefore zs values are fitted to an expression of the form (n + ((constant))^(2/3)) to see if the eigenvalues will conform to the WKB result.   b/100 is used in the fit so that the constant stays in the region where the function is real.

nlm =NonlinearModelFit[Transpose[{Range[0,199],zs}],{a (((n+(b/100)))^(2/3))},{a,b},n]

The fit is quite good.

nlm[{"BestFit","ParameterTable"}]

which results in (Subscript[E, n]) = 1.7708 * (n ((+0.488))^(2/3)) (((ℏ k)/(Sqrt[2m])) )^(2/3), quite close to the WKB result.

The WKB prediction is quite good, even for the ground state.

{zs[[1]], N @ scaledenergies[1]/. n -> 0}

The value of the WKB energy divided by the actual energy goes rapidly to 1 as the quantum number increases.

ListPlot[Thread[{Range[0,9],((Table[scaledenergies[1],{n,0,9}])/(zs[[Range[10]]]))}],PlotRange -> All,PlotStyle -> PointSize[Large], AxesOrigin -> {0,1}, Frame -> True, FrameLabel -> {"n", "WKB/Actual"} ]

Rigid Square Well Potential

A rigid square well potential with walls at -a and a can be regarded as the limit of (((x/a))^p) as p -> ∞

Plot[(((x/2))^100),{x,-3,3}, PlotRange -> {0,5}]

The first step is to replace k by (a^(-p))

energies[p]/. k -> (a^(-p))//PowerExpand

and then take the limit as p -> ∞

Limit[%, p -> ∞]//Simplify

The result obtained(( ((( (n + 1/2) π ℏ))^2))/(8 m (a^2))) is close to the Schrodinger equation result of ((  ((((n+1)  π ℏ))^2))/(8 m (a^2))), where n is the number of nodes (see Shiff, p. 39)

However, the ground state energy (n = 0) obtained from the WKB approximation is 1/4 the actual energy.

The energy scaling factor is

Limit[(((k ((((ℏ/(Sqrt[2 m]))))^p)))^(2/(2+p)))/. k -> (a^(-p))//PowerExpand,p -> ∞]

Therefore the actual ground state energy in dimensionless units is

((((π ℏ))^2)/(8 m (a^2)))/((ℏ^2)/(2 (a^2) m))

N[%]

Numerical Solutions

Quartic Potential

The quartic potential is analyzed first since there are literature results on the ground state energy.

N @ scaledenergies[4]/. n -> 0

The WKB approximation predicts 0.867145 for the scaled energy of the ground state for the quartic potential.  The ground state, being symmetric under inversion about y = 0, has a non-zero value, but zero slope at y = 0.  The value of the wavefunction, as a function of energy, can be found using DSolve, in the form of a DifferentialRoot, which can be evaluated numerically for infinite precision values of e and finite precision values of y.

dr[e_] = ψ /.DSolve[ψ''[y] +( e-(y^4))ψ[y]==0 && ψ[0]==1 && ψ'[0]==0,ψ,y][[1]]

The value of the scaled energy for the ground state can be approximated by setting ψ[10] equal to 0.

sln = FindRoot[dr[Rationalize[e,(10^(-16))]][10.] == 0,{e,0.87}]

This result compares well to the result 1.060362090484182899647046016692663545515208728528977933216245241695
943563044344421126896299134671703510546244358582525580982763829 obtained by Trott

e4 = Rationalize[e /. sln,(10^(-16))]

Cubic Potential

The WKB approximation predicts 0.920791 for the scaled energy of the ground state for the quartic potential.

N @ scaledenergies[3]/. n -> 0

The ground state energy is found using the same process as was used for the quartic potential.

Clear[dr]

dr[e_] = ψ /.DSolve[ψ''[y] +( e-(y^3))ψ[y]==0 && ψ[0]==1 && ψ'[0]==0,ψ,y][[1]]

sln = FindRoot[dr[Rationalize[e,(10^(-16))]][10.] == 0,{e,0.82}]

e3 = e/. Rationalize[sln, (10^(-16))]

Quintic Potential

The WKB approximation predicts 0.8207 for the scaled energy of the ground state for the quartic potential.

N @ scaledenergies[5]/. n -> 0

Clear[dr]

dr[e_] = ψ /.DSolve[ψ''[y] +( e-(y^5))ψ[y]==0 && ψ[0]==1 && ψ'[0]==0,ψ,y][[1]]

sln = FindRoot[dr[Rationalize[e,(10^(-16))]][10.] == 0,{e,0.83}]

e5 = e/. Rationalize[sln, (10^(-16))]

Summary

WKB ground state results in dimensionless units for p = 1 to 5 and ∞

wkbgssc =N @ Join[Table[scaledenergies[p],{p,1,5}], {Limit[scaledenergies[p], p -> ∞]}]/. n -> 0

Actual ground state results in dimensionless units for p = 1 to 5 and ∞

actgssc = N @ {zs[[1]],1 ,e3, e4, e5,2.4674}

The WKB energies divided by the actual energies for p = 1, 2, 3, 4, 5,  ∞ are

(wkbgssc/actgssc)

For powers of 1 to 5, the WKB ground state divided by the actual ground state energy is close to linear in p

ps = {1,2,3,4,5,∞};

ListPlot[Thread[{ps,(wkbgssc/actgssc)}],PlotStyle -> PointSize[Large]]

Inversion of the WKB Quantization Condition

(See Engelke for a slightly different derivation of the inversion.)

If n is regarded as a function of E, the quantization condition

(∫_0^E)(((  E-V  ))^(-(1/2)))x(V)dV =  ( n(E) + (1/2) )(( π ℏ)/(Sqrt[2m]))

is an example of Abel’ s Integral equation.

(∫_0^t)(((  t-τ ))^(-α))u(τ)dτ = f(t) for t > 0

The general solution (Andrews and Schivamoggi, p. 242) is

u(t) =(( sin(π α))/π) (∫_0^t)(((  t-τ ))^(α-1))f'(τ)dτ

Therefore, since α = (1/2)

x(E) = (1/π)(( h)/(2(Sqrt[2m])))(∫_0^E)  n'(V)  (((  E-V  ))^(-(1/2)))dV

Since V is a function of x and n is a function of E, interchanging V and E in the expression gives a result which is more physically meaningful.

x(V) = (1/π)(( h)/(2(Sqrt[2m])))(∫_0^V)  n'(E)  (((  V-E  ))^(-(1/2)))dE

using ℏ instead of h gives

x(V) = (ℏ/(Sqrt[2m]))(∫_0^V)  n'(E)  (((  V-E  ))^(-(1/2)))dE

Infinite Square Well

First find n as a function of E using the analytic result

n[e_] = n/. Solve[(((((n+1))^2)(h^2))/(32 m (a^2)))==e,n][[2]]

and use that to find x(V)

x[v_]= (1/π)(h/(2(Sqrt[2m])))Integrate[n'[e](((v-e))^(-1/2)),{e,0,v}, Assumptions -> {v > 0}]

Since x(0) = 0, x(V) = a is an infinite square well between a and -a.

Symmetric Linear Potential

Assume that an experiment gives values for the lowest 10 energy levels of a linear symmetric potential and the results are fit to E ∝ (n + ((b))^p), motivated by the WKB approximation.
The fit is to n - 1 since the first eigenvalue is for n = 0.

ff = FindFit[zs[[Range[10]]],a (((n-1+(b/10)))^p),{a,b,p},n]

From the fit, n(E) can be calculated

n[e_] = n /. Solve[e ==  (a (((n+b))^p)/. ff),n][[1]]

and then x (V) from that

x[v_]=Integrate[n'[e](((v-e))^(-1/2)),{e,0,v}, Assumptions -> {v > 0}]

x(V) is close to V, the actual dependence.

Conclusions

The WKB quantization condition is remarkably accurate for symmetric power law potentials V(x) = k Abs( ((x ))^p) for p = 1 to 5.  The predicted ground state energy is about 10% too high for p = 1, exact for p = 2, about 10% too low for p = 3, about 20% too low for p = 4 and about 25% too low for p = 5.  In the limit that p -> ∞, the predicted ground state energy is one-quarter the actual value.  For the symmetric linear potential, the error for states above the ground state decreases rapidly with n, being less than 1% for n = 1. Also for the linear symmetric potential, the inverse quantization condition can be used to estimate the potential dependence with good accuracy, using a small number of low lying states.

References

Andrews, L.C. and Shivamooggi, B.K., “Integral Transforms for Engineers and Applied Mathematicians”, MacMillan Publishing, 1988

Engelke, Ray, “An approximate analytic solution for the spectra inverse problem”, Am. J. Phys. 37(4) 1989 p. 339 - 344

Robinett, Richard W. “Quantum Mechanics”  Oxford University Press, 1997

Schiff, Leonard, “Quantum Mechanics”,  McGraw-Hill, 1968

Trott, Michael, “High-Precision Value for the Quartic Anharmonic Oscillator Ground State”, arXiv:quant-ph/0012147