Divergences in perturbation theory

Quantum Mechanics
                       K. Konishi, G. Paffuti

Find a numerical recursive relation for the Hamiltonian eigenvalues of a quartic oscillator.  Compute the first  50 or 60 orders in perturbation theory and apply the Padè resumation technique to approximate the series.

Introduction

The hamiltonian is defined by

H = ((p^2)/2) +((x^2)/2)+ g ((x^4)/2);

We parametrize the solution of Schrödinger equation for eigenstates of H by

-(1/2)ψ''(x) + ((x^2)/2)ψ + g ((x^4)/2)ψ = E ψ;        ψ = B(x) Exp[-((x^2)/2)];

The resulting equation for B[x] is

B'' - 2x B'- g (x^4) B + (2 E - 1 ) B = 0

B can be expanded in g

B[x]= (∑_(k=0)^∞)(g^k)(Subscript[B, k])[x];

By substitution in () it follows that the functions (Subscript[B, k])[x] are polynomials

(Subscript[B, k])[x]=(∑_(j=0)^(2k))(Subscript[A, k,j])(x^(2j))  ; (Subscript[A, k,0])= 1;       2 E = (Underoverscript[∑ , (k=0), ∞])(Subscript[ϵ, k]) (g^k) ;

where coefficients (Subscript[A, k,j]) satisfy

(Subscript[A, k,j])= (1/(4j))((2j+2)(2j+1)(Subscript[A, k,j+1]) - (Subscript[A, k-1,j-2])+(Underoverscript[∑ , (s=1), (k-1)])(Subscript[A, k-s,j]) (Subscript[ϵ, s]));          (Subscript[A, k,2k+1])=0;

This recursion relation can be solved starting with the sequence j = 2k, 2k-1 , ...  1. Once (Subscript[A, k1])is found (Subscript[ϵ, k]) can be computed by

(Subscript[ϵ, k]) = - 2 (Subscript[A, k1 ])- (Underoverscript[∑ , (s=1), (k-1)])(Subscript[ϵ, s]) .

Numerical computation of the perturbation expansion for the ground state

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

The function CofEne just reproduces the above recursion relation. We add a shift in the indexing of the matrix (Subscript[A, kj]) just to avoid negative values in k and j. The variable norder is the requested maximum power of g.

CofEne[norder_]:=Module[{A,Ishift,ene}, 
Ishift = 3; 
(* Ishift è introdotto per comodità *)
A = Table[0,{i,1,norder+2+Ishift},{j,1,2*norder+3+Ishift}];
ene = Table[0,{i,1,norder}];
A[[Range[Ishift,Ishift+norder],Ishift]] = 1;
Do[
Do[
A[[k+Ishift,j+Ishift]] = 1/(4*j)*(
(2j+2)(2j+1)A[[k+Ishift ,j+1+Ishift]] - 
A[[k+Ishift-1,j+Ishift-2]] +
ene[[Range[k-1]]].A[[Range[Ishift + k-1,Ishift+1,-1],j+Ishift]]);,
{j,2*k,1,-1}];
ene[[k]] = - 2*A[[k+Ishift,1+Ishift]] - Total[ene[[Range[1,k-1]]]];,
{k,1,norder}];
(1/2)ene]

We first fix a maximum order for the series

norder =60;

then we compute perturbation coefficients for the energy

ene = CofEne[norder];

maxOrderToBePrinted = 20;
ek=Table[ToString["E"[k],StandardForm]<>" =  "<>ToString[ ene[[k]],StandardForm],{k,1,norder}];
TableForm[Table[{ek[[k]],ek[[k+1]]},{k,1,maxOrderToBePrinted-1,2}]]

We can compare the result with asymptotic predictions (see text) :

eneas = Table[(Sqrt[(6/(π^3))])Gamma[k+(1/2)]((((3/2)))^k),{k,1,norder}];

Show[
 ListPlot[Log[Abs[ene]],PlotRange->{{0,norder+10},{0,Automatic}},
PlotStyle-> PointSize[0.012],AxesLabel->{"n","Log[(Subscript[E, n])]"}],
 ListPlot[Log[Abs[eneas]],
PlotRange->{{0,norder+10},Automatic},PlotStyle-> RGBColor[1,0,0],Joined-> True]]

It is apparent the behavior (Subscript[E, n])∼ n! .  Below we plot the difference between asymptotic prediction and exact values

ListPlot[Log[eneas]-Log[Abs[ene]],AxesLabel->{"n",""},
PlotLabel-> "Log[|(Subscript[E, n])|]: differences"]

Borel-Padé resummation

Assuming some analyticity properties (see text) we can reconstruct a function string from its asymptotic series using a Borel transform:

E[g] = (1/2)+ (1/g)(∫_0^∞)Exp[-(u/g)] B[u]du

From () it is easily seen that order by order in g, E[g] is given by the Taylor expansion of B[u], with

B[u] = (∑_(k=1))(Subscript[E, k])((u^k)/(k!));     E[g] = (∑_(k=1))(Subscript[E, k])(g^k)

For a truncated expansion we have at our disposal only a polynomial approximation of the function B[u]:

B[u] = (∑_(k=1)^(2N))(Subscript[E, k])((u^k)/(k!));

while to perform in a nonperturbative estimation of  () we need the whole B[u]. This function can be obtained by using a Padé approximation, i.e. writing

B[u] = Padè[N,N,u] = ((Polin. di grado N)/(Polin. di grado N))

The following function compute the resummed value of energy for a list gV of values of the coupling constant g:

calcPade[orderPade_,gV_,ene_]:= Module[{borel,borelPade,energy},
borel[z_]:=Sum[ene[[i]]*z^i/i!,{i,1,2*orderPade}];
borelPade= PadeApproximant[borel[z],{z,0,orderPade}];
energy = Table[(1/2)+ (1/(gV[[i]])) NIntegrate[Exp[-z/gV[[i]]]*borelPade,{z,0,∞}],

{i,1,Length[gV]}]]

Let us compute the approximants for different values of the coupling and for different order of the approximants

gValues = {1.,2.,10.};

resummedEnergies = Table[calcPade[ordPade,gValues,ene],{ordPade,10,Floor[norder/2],10}]

We want to compare these values with those obtained by a numerical solution of the Schrödinger equation. To perform the comparison we will use a couple of functions which will be fully explained in the chapter on variational approximation: essentially we take a finite dimensional Hilbert space and we diagonalize exactly the Hamiltonian in that subspace. The reader can consider the function variationalEnergies[g] just as a black-box which computes a variational approximation to eigenvalues.

Variational computation

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

variationalEnergies[gV_]:= Module[{parity,dim,p2,x2,x4,h,α,g,eneTh},
parity = 0; dim = 400;
eneTh = Table[0.,{i,1,Length[gV]}];
{p2,x2,x4} = hermiteMatrices[dim,parity];
Do[ 
g = gV[[i]];
α=x/.FindRoot[(2dim+1)(x-(1./(x^3)))-((3g)/(x^5))(2(dim^2)+2dim+1),{x,(((3g))^(1/3))}];h =((α^2)/2)p2 +(1./(2(α^2) ))x2 + g (1/(2(α^4) ))x4;
eneTh[[i]] = Sort[Eigenvalues[h]][[1]];,
{i,1,Length[gV]}];eneTh]

Comparison of the results

We get variational approximation of the eigenvalue :

eneTh= variationalEnergies[gValues];

Here we print the comparison. It is apparent the agreement.

gString = Table["g="<>ToString[gValues[[i]]],{i,1,Length[gValues]}];
ordString =Append[ToString/@ Range[10,Floor[norder/2],10],"var."];

NumberForm[TableForm[Append[resummedEnergies,eneTh],TableSpacing->{2,4},
TableHeadings-> {ordString,gString}],
{13,12},DigitBlock->3,NumberSeparator->" ",NumberPadding-> {"","0"}]

The agreement improve with a longer perturbation input. The following lines can take some time:

norder =120; 
ene = CofEne[norder];

gValues = {1.,2.,10.};

resummedEnergies = Table[calcPade[ordPade,gValues,ene],{ordPade,10,Floor[norder/2],10}]

eneTh= variationalEnergies[gValues];

Here we print the comparison. It is apparent the improvement

gString = Table["g="<>ToString[gValues[[i]]],{i,1,Length[gValues]}];
ordString =Append[ToString/@ Range[10,Floor[norder/2],10],"var."];

NumberForm[TableForm[Append[resummedEnergies,eneTh],TableSpacing->{2,4},
TableHeadings-> {ordString,gString}],
{13,12},DigitBlock->3,NumberSeparator->" ",NumberPadding-> {"","0"}]