WKB approximation

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study the WKB approximation for anharmonic potential by writing the solution in terms of elliptic integrals.

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

Introduction

This notebook is written to help the reader in the manipulation of elliptic integrals which appear often in the study of anharmonic oscillator. A detailed account of analytical computed integrals for anharmonic oscillator, scattered in several notebooks and problems, is given.

We give also a look at the weak coupling expansion of the results and we compare them with a few orders of perturbation theory.

Anharmonic oscillator μ = 1, g > 0.

Lowest order WKB

The potential is

V[x] = (1/2)(x^2)+ g (1/2)(x^4);

Defining

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

the basic quartic polynomial 2 (E - V) is expressed as

g((a^2)-(x^2))((b^2)+ (x^2))

$Assumptions = {ϵ>0,g>0};
Simplify[g((a^2)-(x^2))((b^2)+ (x^2))== 2(ϵ - V[x])]

x = ± a are the classical turning points

The relevant quantities are the the action integral and the classical period of motion

J[E] = ((2 (Sqrt[g]))/(π ))(∫_0^a)(Sqrt[((a^2)-(x^2))((b^2)+(x^2))]) 𝕕x

T = ∮ ((dx)/(p[x]))= (4/(Sqrt[g]))(∫_0^a)((dx)/(Sqrt[((a^2)-(x^2))((b^2)+(x^2))])).

Elliptic integrals

These integrals, and others considered below can be expressed through elliptic integrals. Their definition is

EllipticE[m] = (∫_0^(π/2))(Sqrt[1- m ((Sin[φ])^2)]) 𝕕φ ;    EllipticK[m]= (∫_0^(π/2)) ((𝕕φ)/(Sqrt[1- m ((Sin[φ])^2)])) .

A transformation property wich will be useful below is the following (all parameters are supposed to give rise to real integrals, as needed in different integrals)

EllipticE[-m] = (Sqrt[1+m])EllipticE[(m/(1+m))];

EllipticK[-m]= (1/(Sqrt[1+m]))EllipticK[(m/(1+m))]

These property are easily shown with the change of variable φ → π/2 - φ.

EllipticE[-m] =(∫_0^(π/2))(Sqrt[1+ m ((Sin[φ])^2)]) 𝕕φ (Underscript[=, φ→π/2-φ]) (∫_0^(π/2))(Sqrt[1+ m ((Cos[φ])^2)]) 𝕕φ =
(∫_0^(π/2))(Sqrt[1+ m - m ((Sin[φ])^2)]) 𝕕φ = (Sqrt[1+m])EllipticE[(m/(1+m))]

and similarly for EllipticK. We write a rule for this transformation

ruleElliptic= {EllipticE[-m_]->  (Sqrt[1+m])EllipticE[(m/(1+m))],
EllipticK[-m_]-> (1/(Sqrt[1+m]))EllipticK[(m/(1+m))]};

Action and period

Clear[a,b];

The action integral is computed as

$Assumptions = {ϵ>0,g>0,a>0,b>0};
Integrate[(Sqrt[((a^2)-(x^2))((b^2)+(x^2))]),{x,0,a}]

j[a_,b_]:= ((2 (Sqrt[g]))/(π ))(1/3) b (((a^2)-(b^2)) EllipticE[-((a^2)/(b^2))]+((a^2)+(b^2)) EllipticK[-((a^2)/(b^2))]);

while the integral for classical period is

Integrate[1/(Sqrt[((a^2)-(x^2))((b^2)+(x^2))]),{x,0,a}]

t[a_,b_]:=(4/(Sqrt[g]))((EllipticK[-((a^2)/(b^2))])/b)

These forms are slightly different from the ones used in notebooks (for hystorical reasons):

((2 (Sqrt[g]))/(π ℏ))((Sqrt[(a^2)+(b^2)])/3)((b^2)EllipticK[((a^2)/((a^2)+(b^2)))] - ((b^2)-(a^2))EllipticE[((a^2)/((a^2)+(b^2)))])
(4/(Sqrt[g]))(1/(Sqrt[(a^2)+(b^2)]))EllipticK[((a^2)/((a^2)+(b^2)))]

They are seen to be equivalent by applying the rule above :

Simplify[j[a,b]/.ruleElliptic]

Simplify[t[a,b]/.ruleElliptic]

Series expansions

In several problems power series in g of the quantities just computed are needed. The most simple way to obtain this series is by using the Mathematica function Series. By reintroducing explicit values for a and b:

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

This is the perturbative series for the action integral and period are

order = 4;

Jofg=Normal[Series[j[a,b],{g,0,order}]]

Tofg = Normal[Series[t[a,b],{g,0,order}]]

Energy quantization

Energy levels are computed by inverting the relation

J = ℏ(n+(1/2))

We write energy as a power series in g:

E = (Subscript[E, 0]) + (∑_k) (g^k) (Subscript[E, k])

The coefficients can be found by inverting the quantization condition.

At leading order (Subscript[E, 0])= ℏ (n+1/2) ≡ s. We write an arbitrary power series

av = Array[c,order];
sub =s+ Sum[(g^k) c[k],{k,1,order}];

Here we extract the coefficient list of the varius power of g once we substitute the series in ℏ J

cf=CoefficientList[s -  Jofg/.ϵ-> sub,g];

The first coefficient is the leading order, independent of g and has already been fixed. In the rest all powers of g must cancel, then

Solve[cf[[2;;order+1]]== 0,av[[1;;order]]]

energy = sub/.First[%]

WKB corrections

ClearAll[a,b];

The next to leading correction to action integral is determined by

J = (1/π)(∫_(Subscript[x, 1])^(Subscript[x, 2]))p[x]dx ;   δJ = - ((ℏ^2)/(24π))(∂/(∂E))(∫_(Subscript[x, 1])^(Subscript[x, 2]))((V''[z])/(((2m(E-V[x])))^(1/2))) .

and the corresponding quantization condition is

J + δJ = n+ (1/2)

Then, in the present case

δJ = - ((ℏ^2)/24)(∂/(∂E))(∫_(-a)^a)((1+ 6g (x^2))/((Sqrt[g])(Sqrt[((a^2)-(x^2))((b^2)+ (x^2))]))) = - ((ℏ^2)/12)(1/(Sqrt[g]))(∂/(∂E))(∫_0^a)((1+ 6g (x^2))/(Sqrt[((a^2)-(x^2))((b^2)+ (x^2))]))

Mathematica can do the integral easily

Integrate[((1+ 6g (x^2))/(Sqrt[((a^2)-(x^2))((b^2)+ (x^2))])),{x,0,a}]

and inserting the definitions of a and b the derivative can be computed. To get a more readeable result we leave a and b in the result and use

(∂_ϵ)a = (1/(a g))(1/((a^2)+(b^2)));  (∂_ϵ)b = (1/(b g))(1/((a^2)+(b^2)));

expr1=FullSimplify[- ((ℏ^2)/(12π))(1/(Sqrt[g]))Dt[(6 b g EllipticE[-((a^2)/(b^2))]+(((1-6 (b^2) g) EllipticK[-((a^2)/(b^2))])/b)),ϵ]/.{Dt[g,ϵ]-> 0,Dt[a,ϵ]-> 1/(a g ((a^2)+(b^2))),Dt[b,ϵ]-> 1/(b g ((a^2)+(b^2)))}]

δj[a_,b_]:= ((ℏ^2) ((-(b^2)+(a^2) (1-12 (b^2) g)) EllipticE[-((a^2)/(b^2))]+((a^2)+(b^2)) (1+6 (a^2) g) EllipticK[-((a^2)/(b^2))]))/(12 (a^2) b ((((a^2)+(b^2)))^2) (g^(3/2)) π)

WKB and perturbation theory

It can be interesting to compute also the perturbative expansion for δJ

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

δJofg=Normal[Series[δj[a,b],{g,0,order}]]

The inversion of the quantization condition is now (we have transferred a factor ℏ to the right hand side to form s = ℏ(n+1/2) )

av = Array[c,order];
sub =s+ Sum[(g^k) c[k],{k,1,order}];
cf=CoefficientList[s - Jofg -  δJofg/.ϵ-> sub,g];
energy1=sub/.First[Solve[cf[[2;;order+1]]== 0,av[[1;;order]]]]

Explicitly, in terms of the principal quantum number

energy1ofg=Collect[Expand[energy1/.s-> ℏ(n+(1/2))],{g,ℏ},Simplify]

to be compared with previous result

energyofg=Collect[Expand[energy/.s-> ℏ(n+(1/2))],{g,ℏ},Simplify]

Both series are different from perturbation theory :

Epert = ℏ(n+(1/2))+ (3/8) g (1+2 n+2 (n^2)) (ℏ^2)-(1/32) (g^2) (21+59 n+51 (n^2)+34 (n^3)) (ℏ^3)+(3/128) (g^3) (111+347 n+472 (n^2)+250 (n^3)+125 (n^4)) (ℏ^4)-(1/2048)(g^4) (30885+111697 n+160470 (n^2)+142610 (n^3)+53445 (n^4)+21378 (n^5)) (ℏ^5);

The k power of g has power of n up to k+1, e.g. (g^3)has powers up to (n^4). For the differences in the above approximations we have

Collect[Simplify[energyofg - Epert],g]

i.e. zeroth order WKB has the first two leading powers in n coincident with the perturbative series. For first order

Collect[Simplify[energy1ofg - Epert],g]

i.e. first order WKB has the first four leading powers in n coincident with the perturbative series.

Anharmonic oscillator μ = -1, g > 0.

ClearAll[a,b,ϵ,g,j,t,δj];

The case of bistable oscillator is a bit more complicated.  The potential is

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

and the 4 inversion points are  x = ± a, ± b where

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

We assume -1/(8g) < E < 0, i.e. 4 real turning points.

The basic quadratic form is

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

In effect

$Assumptions = {1+8ϵ g > 0,g>0};

V[x_]:=-(1/2) (x^2) + (1/2) g (x^4);
a= (Sqrt[((1-(Sqrt[1+8 ϵ g]))/(2g))]);  b= (Sqrt[((1+(Sqrt[1+8 ϵ g]))/(2g))]);

FullSimplify[Expand[2(ϵ - V[x])]==Expand[ g((x^2)-(a^2))((b^2)- (x^2))]]

A plot of the potential is given below :

Block[{ϵ=-1/2,g=0.1},
Plot[V[x],{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]} ]]

Clear[a,b]

In this problem there are two main question, for negative enegy bound states: the semiclassical quantization inside one of the well, and the tunneling effect which mixes the states. As it is well known the tunneling produces a doublet for each level, the energy separation being proportional to the barrier penetration integral, indicated with K below:

Δ E ∼ ℏ (ω/π)K = ℏ (ω/π) Exp[((-1)/ℏ)(∫_a^a)|p[x]|dx

ω = 2 π/T is the classical oscillation frequency inside a well.

In the problems of this chapter the value of the relevant integrals has been given:

J = ((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))]);

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

K = Exp[-2 D];  D = (Sqrt[g])(∫_0^a)(Sqrt[( (a^2)- x^2)((b^2)-(x^2))])dx = (Sqrt[g]) (b/3)(((a^2)+(b^2))EllipticE[((a^2)/(b^2))] -  ((b^2) - (a^2))EllipticK[((a^2)/(b^2))])

Our first task is to derive these results with Mathematica .

Mathematica refuses to perform the integrals in this form

TimeConstrained[Integrate[(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]),{x,a,b},Assumptions-> {0<a<b}],3]

A possible help is to write

(b^2)= (a^2)+ (c^2);   (x^2) = (a^2)+ (c^2) (z^2);

as z varies from 0 and 1 x span the interval (a,b). The jacobian is given by

x dx = (c^2)zdz ⇒ dx = ((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)]))z dz

and we have

(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = c z (Sqrt[((b^2)- (a^2)- (c^2)(z^2))]) ((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)])) z dz =((c^4)/c)((Sqrt[1- (z^2) ])/(Sqrt[((a^2)/(c^2))+ (z^2)]))(z^2)dz

Now Mathematica  perform the integral

((c^3)Integrate[((Sqrt[1- (z^2) ])/(Sqrt[((a^2)/(c^2))+ (z^2)]))(z^2),{z,0,1},Assumptions-> {a>0,c>0}])/.{(c^2)-> (b^2)- (a^2)}

j[a_,b_]:= ((Sqrt[g])/π)(1/3) a (((a^2)+(b^2)) EllipticE[-((-(a^2)+(b^2))/(a^2))]-2 (b^2) EllipticK[-((-(a^2)+(b^2))/(a^2))])

This integral can be reported to the form printed above with the rule written at the beginning of the notebook :

Simplify[j[a,b]/.ruleElliptic,Assumptions-> {a>0,b>0}]

The same trick can be used to compute the classical period

(1/(Sqrt[((x^2)- (a^2))((b^2)-(x^2))]))dx = (1/(c z (Sqrt[((b^2)- (a^2)- (c^2)(z^2))]))) ((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)])) z dz =
=(1/c)(1/((Sqrt[1- (z^2) ])(Sqrt[((a^2)/(c^2))+ (z^2)])))dz

((1/c)Integrate[(1/((Sqrt[1- (z^2) ])(Sqrt[((a^2)/(c^2))+ (z^2)]))),{z,0,1},Assumptions-> {a>0,c>0}])/.{(c^2)-> (b^2)- (a^2)}

t[a_,b_]:=(2/(Sqrt[g])) ((EllipticK[-((-(a^2)+(b^2))/(a^2))])/a)

Again ruleElliptic shows the equivalence with the quoted result :

Simplify[t[a,b]/.ruleElliptic,Assumptions-> {a>0,b>0}]

For the transmission integral D there are no particular problems

Integrate[(Sqrt[( (a^2)- x^2)((b^2)-(x^2))]),{x,0,a},Assumptions-> {0<a<b}]

and the expression is the one used in problems and notebooks.

d[a_,b_]:=  (Sqrt[g])(1/3) b (((a^2)+(b^2)) EllipticE[((a^2)/(b^2))]+((a^2)-(b^2)) EllipticK[((a^2)/(b^2))])

First order WKB

The correction to the action integral

δJ = - ((ℏ^2)/(24π))(∂/(∂E))(∫_(Subscript[x, 1])^(Subscript[x, 2]))((V''[z])/(((2m(E-V[x])))^(1/2))) .

takes the form (we put ℏ = 1)

δJ = - (1/(24π))(1/(Sqrt[g]))(∂/(∂E))(∫_a^b)((-1 + 6 g (x^2))/(Sqrt[((x^2)-(a^2))((b^2)- (x^2))])) .

Let us perform the change of variables as above

(Subscript[((-1 + 6 g (x^2))/(Sqrt[((x^2)-(a^2))((b^2)- (x^2))]))dx = ((-1 + 6 g ( (a^2)+ (c^2) (z^2)))/(c z c(Sqrt[1-(z^2)])))((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)]))z dz = ((-1 + 6 g ( (a^2)+ (c^2) (z^2)))/(c(Sqrt[1-(z^2)])(Sqrt[((a^2)/(c^2))+  (z^2)])))dz =
((-1 + 6 g ( (a^2)+ (c^2) ))/(c(Sqrt[1-(z^2)])(Sqrt[((a^2)/(c^2))+  (z^2)])))dz- 6g c  ((Sqrt[1-(z^2)])/(Sqrt[((a^2)/(c^2))+  (z^2)]))= ((-1 + 6 g (b^2))/(c(Sqrt[1-(z^2)])(Sqrt[((a^2)/(c^2))+  (z^2)])))dz - 6g c  ((Sqrt[1-(z^2)])/(Sqrt[((a^2)/(c^2))+  (z^2)]))dz=
((-1 + 6 g (b^2))/(c(Sqrt[1-(z^2)])(Sqrt[((a^2)/(c^2))+  (z^2)])))dz - - 12g c (∂_β)(Sqrt[1-(z^2)])(Sqrt[β+  (z^2)])|, β = ((a^2)/(c^2))])

Let us perform separately the two integrals

Clear[a,b]

int1=Integrate[((-1 + 6 g (b^2))/(c(Sqrt[1-(z^2)])(Sqrt[((a^2)/(c^2))+  (z^2)]))),{z,0,1},Assumptions->{a>0,c>0}]

int2= Simplify[-12g c (∂_β)Integrate[(Sqrt[1-(z^2)])(Sqrt[β+  (z^2)]),{z,0,1},Assumptions->{β>0}]/.β-> ((a^2)/(c^2)),Assumptions-> {a>0,c>0}]

and finally

totInt = Simplify[(int1+int2)/.{(c^2)-> (b^2)- (a^2)}]

Now we have to perform the derivative with respect to energy. We reintroduce the values of a an b using the derivatives:

(∂_ϵ)a =- (1/(a g))(1/((b^2)-(a^2)));  (∂_ϵ)b = (1/(b g))(1/((b^2)-(a^2)));

Clear[a,b];

expr1=FullSimplify[- (1/(24π))(1/(Sqrt[g]))Dt[(6 a g EllipticE[1-((b^2)/(a^2))]-((EllipticK[1-((b^2)/(a^2))])/a)),ϵ]/.{Dt[g,ϵ]-> 0,Dt[a,ϵ]-> -1/(a g (-(a^2)+(b^2))),Dt[b,ϵ]-> 1/(b g (-(a^2)+(b^2)))}]

δj[a_,b_]:= (((b^2)+(a^2) (1-12 (b^2) g)) EllipticE[1-((b^2)/(a^2))]+2 (b^2) (-1+3 ((a^2)+(b^2)) g) EllipticK[1-((b^2)/(a^2))])/(24 a (b^2) ((((a^2)-(b^2)))^2) (g^(3/2)) π);

The reader can verify that the same result would be obtained by writing explicitly

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

Perturbative expansion

Leading order WKB

It is convenient to shift energies, putting the zero at the bottom of the potential

ϵ → -1/(8 g) + ϵ

In this way :

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

The perturbative expansion of various quantities is easily computed :

order = 2;

Jofg=Normal[Series[j[a,b],{g,0,order}]]

δJofg=Simplify[Normal[Series[δj[a,b],{g,0,order}]]]

Tofg = Simplify[Normal[Series[t[a,b],{g,0,order}]]]

Dofg = Expand[Simplify[Normal[Series[d[a,b],{g,0,2}]],Assumptions-> {g>0,ϵ>0}]]

Anharmonic oscillator μ = 1, g < 0.

Clear[a,b,ϵ,V]

This is the case of a metastable system, let us write the potential in the form

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

The classical turning points are ± a and ± b where

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

The basic quartic polynomial is factorized in terms of a and b as

(p^2)[x] = 2m(ϵ - V[x]) = g((a^2)-(x^2))((b^2)- (x^2))

$Assumptions = {ϵ>0,g>0};
Simplify[g((a^2)-(x^2))((b^2)- (x^2))== 2(ϵ - V[x])]

Clear[a,b]

In the problems of this chapter the value of the relevant integrals has been given:

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))])

The first two integral are easily computed

Clear[j,d,t,δj]
$Assumptions = {ϵ>0,g>0,a>0,b>0};
(2((Sqrt[g])/π))Integrate[(Sqrt[((a^2)-(x^2))((b^2)-(x^2))]),{x,0,a},Assumptions-> 0<a<b]

j[a_,b_]:= ((2 b (Sqrt[g]) (((a^2)+(b^2)) EllipticE[((a^2)/(b^2))]+((a^2)-(b^2)) EllipticK[((a^2)/(b^2))]))/(3 π))

(2/(Sqrt[g]))Integrate[(1/(Sqrt[((a^2)-(x^2))((b^2)-(x^2))])),{x,-a,a},Assumptions-> 0<a<b]

t[a_,b_]:= ((4 EllipticK[((a^2)/(b^2))])/(b (Sqrt[g])))

Mathematica needs some help for third integral. With

(b^2)= (a^2)+ (c^2);   (x^2) = (a^2)+ (c^2) (z^2);

as z varies from 0 and 1 x span the interval (a,b). The jacobian is given by

x dx = (c^2)zdz ⇒ dx = ((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)]))z dz

and we have

(Sqrt[((x^2)- (a^2))((b^2)-(x^2))])dx = c z (Sqrt[((b^2)- (a^2)- (c^2)(z^2))]) ((c^2)/(Sqrt[(a^2)+ (c^2) (z^2)])) z dz =((c^4)/c)((Sqrt[1- (z^2) ])/(Sqrt[((a^2)/(c^2))+ (z^2)]))(z^2)dz

We have in effect already computed this integral :

risD=((c^3)Integrate[((Sqrt[1- (z^2) ])/(Sqrt[((a^2)/(c^2))+ (z^2)]))(z^2),{z,0,1},Assumptions-> {a>0,c>0}])/.{(c^2)-> (b^2)- (a^2)}

which is equivalent to the one quoted above

Simplify[risD/.ruleElliptic,Assumptions-> 0<a<b]

d[a_,b_]:= (Sqrt[g])(1/3) b (((a^2)+(b^2)) EllipticE[1-((a^2)/(b^2))]-2 (a^2) EllipticK[1-((a^2)/(b^2))])

For completeness we write also the first order WKB correction to J

δJ = - ((ℏ^2)/(24π))(∂/(∂E))(∫_(-a)^a)((V''[z])/(((2m(E-V[x])))^(1/2))) = - ((ℏ^2)/(12π))(1/(Sqrt[g]))(∂/(∂E))(∫_0^a)((1- 6g (x^2))/(Sqrt[((a^2)-(x^2))((b^2)- (x^2))]))

The integral is given by

Integrate[((1-6g (x^2))/(Sqrt[((a^2)-(x^2))((b^2)-(x^2))])),{x,0,a},Assumptions-> 0<a<b]

Now we have to perform the derivative with respect to energy. We reintroduce the values of a an b using the derivatives:

(∂_ϵ)a = (1/(a g))(1/((b^2)-(a^2)));  (∂_ϵ)b = - (1/(b g))(1/((b^2)-(a^2)));

expr1=FullSimplify[- (1/(12π))(1/(Sqrt[g]))Dt[(6 b g EllipticE[((a^2)/(b^2))]+(((1-6 (b^2) g) EllipticK[((a^2)/(b^2))])/b)),ϵ]/.{Dt[g,ϵ]-> 0,Dt[a,ϵ]-> 1/(a g (-(a^2)+(b^2))),Dt[b,ϵ]->- 1/(b g (-(a^2)+(b^2)))}]

Clear[δj]

δj[a_,b_]:= (((-(b^2)+(a^2) (-1+12 (b^2) g)) EllipticE[((a^2)/(b^2))]+((a^2)-(b^2)) (-1+6 (a^2) g) EllipticK[((a^2)/(b^2))])/(12 b ((((a^3)-a (b^2)))^2) (g^(3/2)) π))

The reader can check that the same result is obtained by substituting

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

and taking the derivative with respect to ϵ.

Perturbative expansions

The perturbative expansions are generated by substituting explicit values for a and b

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

This is the perturbative series for the action integral and period are

order = 2;

Jofg=Normal[Series[j[a,b],{g,0,order}]]

Tofg = Normal[Series[t[a,b],{g,0,order}]]

dofg = Normal[Series[d[a,b],{g,0,0}]]

δJofg=Normal[Series[δj[a,b],{g,0,order}]]