This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Recently a divergence-free generalization of the classical WKB approximation was suggested [✶]. Starting from the classical Schrödinger equation -ψ''(x)+V(x)ψ(x)=ε ψ(x), one obtains after making the ansatz ψ(x)~exp(φ(x)) the Riccati equation φ''(x)+((φ'(x))^2)=(Subscript[𝒦, ε])(x), where (Subscript[𝒦, ε])(x)=ε-V(x). Differentiating this equation repeatedly and substituting the second derivatives back after each differentiation step yields an equation of the form (φ^((N)))(x)+(Subscript[𝓅, N])(φ'(x),(Subscript[𝒦, ε]^((k)))(x)).
The quantization condition of the Nth order generalized WKB approximation is then

(1/(2π i))(∮_𝒞)y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x))dx=n+1-(1/2)(Subscript[δ, N,2])

where y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x)) is a solution of (Subscript[𝓅, N])(y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x)),(Subscript[𝒦, (Subscript[ε, n])]^((k)))(x))=0 and the N-dependence of (Subscript[ε, n]) is not shown explicitly.

The integration path is enclosing two branch points along the real axis, and no other singularities of any kind and the sheet has to be chosen in such a way, that y=y(x) is real-valued as x->∞.

For the harmonic oscillator (Subscript[𝒦, ε])(x)=ε-(x^2), check if the first eight quantization conditions result in the well-known ground state eigenvalue (Subscript[ε, 0])=1.

(* Riccati for of the Schrödinger equation *)
eqSR = φ''[x] + φ'[x]^2 + 𝒦[x];
(* replace second derivative with first and 𝒦 *)
rule = Solve[eqSR == 0, Derivative[2][φ][x]][[1]];
(* n-th order equation *)
eq[𝒩_] := Nest[(D[#, x] /. rule)&, eqSR, 𝒩 - 2] // Expand

Here are the first few differential equations.

Table[eq[𝒩], {𝒩, 2, 6}]

The integrand of the quantization condition is easily obtained from the last equations by ignoring the highest derivative and replacing the first derivative by a new variable.

ae[𝒩_] := Collect[(#/Coefficient[#, y[x], 𝒩])&[
                  eq[𝒩]  /. {Derivative[𝒩][φ][x] -> 0, 
                  φ'[x] -> y[x]}], y[x], Factor]

Here are the first few integrands.

Table[ae[𝒩], {𝒩, 2, 6}] /. y𝒦_[x] -> y𝒦

We now specialize to the harmonic oscillator potential. This is value of 𝒦 for the harmonic oscillator expressed as a pure function.

𝒦Function[ε_, x_] := ε - x^2

These are the resulting harmonic oscillator quantization integrands.

ahoQC[𝒩_, ε_, {ξ_, η_}] := ae[𝒩] /. 
          𝒦 -> Function[x, Evaluate[𝒦Function[ε, x]]] /. {x -> ξ, y[x] -> η}

Table[ahoQC[𝒩, ε, {x, y}], {𝒩, 2, 6}]

For a given ε, we must calculate the branch points and other singularities in the complex x-plane for a given quantization condition integrand of order N (which is just an algebraic curve). Two branch points are on the real axis and are needed for the contour integration. But we also need the locations of the other singularities to make sure that we do not enclose any of them in the contour. The function branchPointsAndPoles calculates the branchpoints and other singularities for a given polynomial in w and z.

branchPointsAndPoles[poly_, {w_, z_}] := 
Module[{bpSol, branchPoints, poly1, singularityCondition, singularities},
 (* branch points as multiple roots that connect sheets *)
 bpSol = Solve[Resultant[poly, D[poly, w], w] == 0, z];
 branchPoints = N[#, 50]& @ If[bpSol == {}, {}, z /. bpSol];
 (* all singularities *)
 poly1 = poly /. w -> w[z];
 singularityCondition = Numerator[Together[(1/D[w[z], z] /. 
                          Solve[D[poly1, z] == 0, D[w[z], z]][[1]])]];
 sSol = Solve[Resultant[singularityCondition, poly1, w[z]] == 0, z];
 singularities = N[#, 50]& @ If[sSol == {}, {}, z /. sSol];
 (* poles are non-branch point singularities *)
 {N[branchPoints, 30], 
  Complement[N[singularities, 30], N[branchPoints, 30]]}]

A quick look at the branch points in the Re(y),Im(y),ε-space. We see multiple branch points along the real axis for larger n.

Show[GraphicsArray[#]]& /@ Partition[
Table[(* complex x-planes as a function of ε vertically *)
Graphics3D[Table[{Hue[0.8 ε/2], Point[{Re[#], Im[#], ε}]& /@ 
                  branchPointsAndPoles[ahoQC[𝒩, ε, {x, y}], {y, x}][[1]]},
                  {ε, 0, 2, 1/50}], PlotRange -> All, BoxRatios -> {1, 1, 1},
          PlotLabel -> " 𝒩 = " <> ToString[𝒩]], {𝒩, 2, 7}], 2]

To calculate the contour integral (Subscript[ℐ, n])=(∮_𝒞)(Subscript[𝓅, n])(y,(Subscript[𝒦, ε]^((k)))(x))dx efficiently, and to avoid any problems with crossing branch cuts, we will not use NIntegrate, but instead solve the differential equation ∂(Subscript[ℐ, n])(t)∂t=y((Subscript[𝒦, ε]^((k)))(x)) along the contour 𝒞 parametrized by t. Because the integrand is a solution of an algebraic equation (with x as the second variable and ε as a parameter), we calculate y((Subscript[𝒦, ε]^((k)))(x)) by solving the differential equation ∂(Subscript[𝓅, N])(y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x)),(Subscript[𝒦, (Subscript[ε, n])]^((k)))(x))/∂x=0 along the curve 𝒞. The two differential equations can be solved simultaneously. This ensures automatically that we are always on the same sheet. No branch cuts even exist in this approach, just branch points which have to be included or excluded in the contour.

Now, we have to carry out integration numerically along an ellipse enclosing (the) two branch points on the real axis and having half axes a and b. We start on the real axis with real argument for y. There might be more than one choice for the y values, as well for the branch points, so we try all possible ones. (The condition (Subscript[y, n])(x->∞)->∞ does not uniquely sort out one sheet from the Riemann surface; this condition can be fulfilled on more than one sheet.) The function makeNDSolveInput sets up the differential equations to be solved.

(* select positive real solutions *)
positivSolutions[poly_, y_] :=  
Select[y /. Solve[poly == 0, y], (Im[#] == 0 && Re[#] > 0)&]

makeNDSolveInput[𝒩_, ε_, {t_, 𝓎_, ℐ_}] :=
Module[{(*branchPoints, poles, realLineBranchPoints, 
        smallestVerticalDistanceSingularity, bE, aEs, ξ, 𝓎Ints, a, b*)},
 (* calculate branch points and potential poles *)
 {branchPoints, poles} = branchPointsAndPoles[ahoQC[𝒩, ε, {x, y}], {y, x}];
  realLineBranchPoints = Select[branchPoints, (Im[#] == 0 && Re[#] > 0)&]; 
  smallestVerticalDistanceSingularity = 
          Min[Abs[Select[Flatten[{branchPoints, poles}], Im[#] != 0&]]];
 (* ellipse parameters so that only the two branch points
    on the real axis are enclosed *)
 bE =  If[smallestVerticalDistanceSingularity == Infinity, 1, 
          smallestVerticalDistanceSingularity/10];
 aEs = Plus @@@ Partition[Append[realLineBranchPoints, 
                                 2 Max[realLineBranchPoints]], 2, 1]/2;
 ξ[τ_] := a Cos[τ] + I b Sin[τ];
 (* differential equations for y_N and quantization integral
    along the ellipse *)
 odes = {D[ℐ[x], x] == 𝓎[x], D[ahoQC[𝒩, ε, {x, 𝓎[x]}], x] == 0} /. 
         Derivative[1][f_][x] :> 1/D[ξ[t], t] Derivative[1][f][t] /.
         {𝓎[x] -> 𝓎[t], ℐ[x] -> ℐ[t]} /. x -> ξ[t];
 (* add initial conditions so that y_N is real positive on the real axis *)
 NDSolveFirstArguments = Flatten[#, 1]& @ 
 (* all possible ellipses *)
 Table[𝓎Ints = positivSolutions[ahoQC[𝒩, ε, {aEs[[k]], y}], y];
       Table[print["{aE, y0} = ", {aEs[[k]], 𝓎Ints[[i]]} // N];
             Join[odes /. {a -> aEs[[k]], b -> bE}, 
                  {ℐ[0] == 0, 𝓎[0] == 𝓎Ints[[i]]}], {i, Length[𝓎Ints]}],
       {k, Length[aEs]}]]

Here is one of the resulting ODE systems for N=3.

makeNDSolveInput[3, 1, {t, y1, I1}][[-1]] // N

It remains to carry out the numerical integration of the differential equations. For a given quantization order, and a given ε, we calculate the quantization integrals. We use high-precision arithmetic to solve the differential equations to ensure the correctness of the result. Calculating the left-hand side of the quantization conditions should, for the ground state (Subscript[ε, 0])=1, yield the result 1/2 for N=2, and 1 else. (If we would not know the energies (Subscript[ε, n]), then we would have to carry out an additional root finding process.)

contourIntegralValues[𝒩_, ε_] :=
Module[{odes, ndsols, y1, I1},
 (* y_N and quantization integral ODEs *)
 odes = makeNDSolveInput[𝒩, ε, {t, y1, I1}];
 (* make sure numerical solution is reliable *)
 ndsols = NDSolve[#, {y1, I1}, {t, 0, -2Pi},
                  PrecisionGoal -> 15, WorkingPrecision -> 25]& /@ odes;
 I1[-2Pi]/(2Pi I) /. ndsols]

Here is the classical WKB case—we obtain the value 1/2.

contourIntegralValues[2, 1] // Chop[#, 10^-8]&

We obtain the value 1 for exactly one pair of branch points for the cases N=3,4,5,6.

contourIntegralValues[3, 1] // Chop[#, 10^-8]&

contourIntegralValues[4, 1] // Chop[#, 10^-8]&

contourIntegralValues[5, 1] // Chop[#, 10^-8]&

contourIntegralValues[6, 1] // Chop[#, 10^-8]&

But the cases N=7,8 do not give a 1 for any of the contours. The “best” contours give a result that is off by a few percent.

contourIntegralValues[7, 1] // Chop[#, 10^-8]&

contourIntegralValues[8, 1] // Chop[#, 10^-8]&