anharmonic oscillator - 2

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the effect of an external linear potential (electric field) on a double well system. Verify the Feynman-Hellmann theorem.

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

Introduction

Functions

We need  only one function, which is just a simplification of a matrix used in notebook [9.1]. The calling sequence is

{p2, x2, x4} = hermiteMatrices[dim, parity]

It produces matrix elements for operators p^2 , x^2, x^4. Parity is 0 or 1 for even and odd states. dim is the dimension of the basis. If parity takes any other value the full matrices are computed, for both even and odd states, and the output is changed adding matrix elements of x in the last place. The default is parity = -1, i.e. full matrices

{p2, x2, x4,x1} = hermiteMatrices[dim, -1],  or  hermiteMatrices[dim].

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[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

This function give fixes the scale in the variational approach

getscale[n_,g_,μ_]:= x/.FindRoot[(2n+1)(x-(μ/(x^3)))-((3g)/(x^5))(2(n^2)+2n+1),{x,(((3g))^(1/3))}]

The base functions are given in general by

fBase[k_,x_]:= (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,x]Exp[-x^2/2]

For fixed parity we can use:

fBase[n_,x_,par_]:=(1/(π^(1/4)))(1/(Sqrt[(2^(2n-2+par))(2n-2+par)!]))HermiteH[2n-2+par,x]Exp[-x^2/2]

Wave function computation

Here we just write the recursion equation of rHermite polynomials. It is compiled for speed and the integer n represents the number of polynomials, i.e. all polynomials from degree 0 to degree n-1.

recursiveHermitePolynomials= Compile[{{n,_Integer},{x,_Real,1}},
Module[{y1,y0,p1,ris,ris0},
y0 = 1.0 + 0.*x;
y1 = 2.x;
ris0={y0,y1};
ris=Table[
p1=2. x y1 - 2.  k y0;
y0=y1;
y1=p1;
p1,{k,1,n-2}];
ris=Join[ris0,ris];
ris]];

This function computes normalized Hermite polynomials, in the form needed for the oscillator wave functions. We write two versions of it , one with arguments [dim, r] gives all polynomials up to degree dim-1.  The second with three arguments, [dim, r, par], gives the first dim polynomials of parity fixed by par, (0 for even polynomials, 1 for odd polynomials).

hermiteP[dim_,xIn_]:= Module[{cf,x,ris,coeff},
x = Flatten[{xIn}];
 coeff=Table[(1/((N[π])^(1/4))) (1/(Sqrt[(2.^(k-1))(k-1)!])),{k,1,dim}]; 
ris= recursiveHermitePolynomials[dim,x]*coeff;
If[Length[xIn]== 0,Flatten[ris],ris]];

hermiteP[dim_,xIn_,par_]:= Module[{cf,dimUp,index,x,ris},
x = Flatten[{xIn}];
dimUp=2 dim - 2 + par+1;
index=Range[par+1,2 dim-2+par+1,2];
ris=hermiteP[dimUp,x][[index]];
ris];

A more efficient version would use the recursion relation involving only even or odd polynomials.  We leave it to as an exercise to write  such a routine.

Finally getPsi gives the wave function for a coefficient list eigv at points r. Parity can be specified, if not all polynomials are considered. This is the default.

getPsi[eigv_,r_,par_: -1]:=
If[par== 1|| par==0, 
(eigv.hermiteP[Length[eigv],r,par])*Exp[-(r^2)/2],
(eigv.hermiteP[Length[eigv],r])*Exp[-(r^2)/2]];

Stark effect

The model

μ=-1; g = 0.1;

Plot[μ (1/2)(x^2) +(1/2) g (x^4) - 0.175 x,{x,-1.5(Sqrt[((Abs[μ])/g)]),1.5(Sqrt[((Abs[μ])/g)])},AxesLabel-> {"x", "V(x)"}]

Energy in external field

dim=50;

α=getscale[dim,g,μ];

{p2,x2,x4,x1} =N/@ hermiteMatrices[dim];

ℰ = Join[Range[0,0.02,0.001],Range[0.025,0.1,0.005]];

Ham0= (1/2)(α^2)p2 + (μ/(2(α^2)))x2 + (g/(2(α^4)))x4;

eneAndDipole=
Block[{HamE,eigH,eigV,ord,dip},
Table[
HamE = N[Ham0 - ℰ[[j]]  x1/α];
{eigH,eigV} = Eigensystem[HamE];
ord = Ordering[eigH];
eigH = eigH[[ord]]; 
eigV = eigV[[ord]];
dip1 =  eigV[[1]]. (x1.eigV[[1]])/α;
dip2 =  eigV[[2]]. (x1.eigV[[2]])/α;
{eigH[[1]],eigH[[2]],dip1,dip2},{j,1,Length[ℰ]}]];

eneAndDipole[[All,1]]

ΔE = Table[eneAndDipole[[All,k]]-eneAndDipole[[1,k]],
{k,1,2}];

dipole = Table[eneAndDipole[[All,k+2]],
{k,1,2}];

lev=1;

ListPlot[Transpose[{ℰ,ΔE[[lev]]}],AxesLabel-> {"ℰ","E(ℰ) - E(0)"}]

fΔE = Interpolation[Transpose[{ℰ,ΔE[[lev]]}]];

Show[
ListPlot[Transpose[{ℰ,dipole[[lev]]}],PlotStyle-> PointSize[0.015],PlotRange-> All],
 Plot[-Derivative[1][fΔE][x],{x,0,0.1},
PlotStyle-> RGBColor[1,0,0],PlotRange-> All],AxesLabel-> {"ℰ","d(ℰ)"}]

Excited state

lev=2;

ListPlot[Transpose[{ℰ,ΔE[[lev]]}],AxesLabel-> {"ℰ","E(ℰ) - E(0)"}]

fΔE = Interpolation[Transpose[{ℰ,ΔE[[lev]]}]];

Show[
ListPlot[Transpose[{ℰ,dipole[[lev]]}],PlotStyle-> PointSize[0.015],PlotRange-> All],
 Plot[-Derivative[1][fΔE][x],{x,0,0.1},
PlotStyle-> RGBColor[1,0,0],PlotRange-> All],AxesLabel-> {"ℰ","d(ℰ)"}]