Unstable anharmonic oscillator

Quantum Mechanics
                       K. Konishi, G. Paffuti

Use the WKB approximation to study the presence of metastable states in the potential

V[x]= (1/2)(x^2)- (1/2)g (x^4) .

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

Introduction

In this notebook we will study the occurrence of metastable states in an instable version of the anharmonic oscillator

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

This model has classical turning points at x= ± a, x = ± b for energy ϵ < 1/8g

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

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

WKB approximation

Analytical results

Using a and b as above we can write, in the classical allowed region

2(E - V[x]) = g((a^2)-(x^2))((b^2)- (x^2))

The WKB quantization condition is

J = ((Sqrt[g])/π)(∫_(-a)^a)(Sqrt[((a^2)-(x^2))((b^2)- (x^2))])dx = n+(1/2)

This equation fix energy as a function of g for a gien quantum number.  The maximum possible energy compatible with a metastable state is E = 1/8g, at which a=b = 1/(Sqrt[2g]). In this limiting case

(Subscript[J, lim]) = ((Sqrt[g])/π)(∫_(-a)^a)(Sqrt[((a^2)-(x^2))((b^2)- (x^2))])dx = ((Sqrt[g])/π)(∫_(-a)^a)((a^2)-(x^2))dx = ((Sqrt[2])/(3 g π)) ∼ (0.15/g)

As this value must be greater than 1/2 to allow a WKB bound state we have the condition g < 0.3.

If g is small enough there exist semiclassical metastable states which decay by tunnel effect through the barrier.

As an order of magnitude the width of the state is

Γ =2(ℏ/T) Exp[-2 D];

where T is the classical period of the bound orbit and K = Exp[-2 D] the penetration barrier factor :

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

The integral involved in these quantities can be expressed through elliptic integrals EllipticE and EllipticK, see notebook [*]

J = ((Sqrt[g])/π)(∫_(-a)^a)(Sqrt[((a^2)-(x^2))((b^2)- (x^2))])dx = ((2(Sqrt[g]))/π)(b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] - ((b^2)- (a^2))EllipticK[((a^2)/(b^2))]);
T = (2/(Sqrt[g]))(∫_(-a)^a)(1/(Sqrt[((a^2)-(x^2))((b^2)- (x^2))]))dx = (4/(Sqrt[g]))(1/b)EllipticK[((a^2)/(b^2))];
Γ =2(ℏ/T) Exp[-2 D];  D = (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))])

A more refined analysis can be made at small coupling, see problem [*], which allow a precise determination of the prefactor in ()

Γ (Underscript[→, g → 0]) (4/(Sqrt[π g]))Exp[((-2)/(3g))]

Numerical results

Clear[a,b,g,ϵ];

As in other notebooks devoted to anharmonic oscillator we will use above analytic expression and at the same time a purely numerical approach, the comparison will provide a check for the analytical results.

Using the same notation of notebook [*] we can write

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

fTx4[ene_,g_]:= Module[{a,b},
a= (Sqrt[((1-(Sqrt[1-8 ene g]))/(2g))]);  b= (Sqrt[((1+(Sqrt[1-8 ene g]))/(2g))]);
 (4/(Sqrt[g]))(1/b)EllipticK[((a^2)/(b^2))]];

fKx4[ene_,g_]:= Module[{a,b,dF},
a= (Sqrt[((1-(Sqrt[1-8 ene g]))/(2g))]);  b= (Sqrt[((1+(Sqrt[1-8 ene g]))/(2g))]);
 dF=(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))]);
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].

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

calcK[ene_,g_]:= Module[{aL,aR,Jval,dF},
aL= (Sqrt[(( 1 - (Sqrt[1- 8 g ene]))/(2g))]);aR=(Sqrt[(( 1 + (Sqrt[1- 8 g ene]))/(2g))]);
dF =Re[NIntegrate[(Sqrt[2(-ene +((1/2)(x^2)- (g/2)(x^4)))]),{x,aL,aR}]];
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,aL}]]];

We prepare a bunch of random couples (E,g) to perform the check. E is chosen between 0 and the maximum 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)]

WKB computation

These functions perform WKB quantization for the model

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

Here is an example

g=0.05;

eneWKB0 = getEneWKBx4[0,g]

ΓWKB0 = (2/(fTx4[eneWKB0,g]))fKx4[eneWKB0,g]

Clear[a,b,g,ϵ]

Small g

At small g the elliptic integrals can be expanded in a Taylor series, with the result (see notebook [*] )

Jofg = ϵ+((3 g (ϵ^2))/4);
Tofg = 2 π+3 g π ϵ;
Dofg = (1/(3 g))+(1/2) (-ϵ-4 ϵ Log[2]+ϵ Log[g]+ϵ Log[2 ϵ]);

Then the "naive" semiclassical expression for Γ has an expansion

Simplify[Normal@Series[((2/Tofg)Exp[-2 Dofg]/.ϵ-> 1/2),{g,0,0}],g>0]

which differs from exact perturbative result for a factor (Sqrt[π/e]), this effect has been analyzed in the case of a double well.

GammaPert[g_]:= (4/(Sqrt[π g]))Exp[((-2)/(3g))]

As an example

listG = Range[0.01,0.1,0.01];

eneWKB0 = getEneWKBx4[0,#1]&/@listG;

ΓWKB0 =( (2/(fTx4@@#1))fKx4@@#1)&/@Transpose[{eneWKB0,listG}];

TableForm[Transpose[{GammaPert[listG],ΓWKB0,(Sqrt[(π/E)])ΓWKB0}]]

This model will be analyzed in more detail using complex scaling techniques in the next chapter.