Spectral density in d=1

Quantum Mechanics
K. Konishi, G. Paffuti

Compute the semiclassical approximation for the spectral density in a harmonic oscillator, a purely quartic oscillator and an anharmonic oscillator.

Compare it with the numerical results.

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

Introduction

Summary of the theory

We will consider in this notebook only system with a purely discreet spectrum.

In problems with a single quantum number n there is a definite correspondence between energy and n, n = n[E]. For these systems the spectral density can be written as

g[E] =(∑_j) δ[E-(Subscript[E, j])]= d[n] (dn/dE)(1 + 2 (∑_(k=1)^∞)Cos[2π k n(E)]).

d[n] is the degeneracy. For one dimensional systems d[n] = 1.

In semiclassical approximation n is related to the action integral I. Each motion is such systems is periodic and from WKB approximation it is known that

(dn/dE) = (1/(ℏ ω)).

With this information the semiclassical approximation to () reads

(Subscript[g, SC])(E) = (Subscript[g, TF])(E)+ (2/(ℏω)) (∑_(k=1)^∞)Cos[2π k ((I/ℏ) - (μ/4))]≡(1/(ℏω))+(2/(ℏω)) (∑_(k=1)^∞)Cos[2π k ((I/ℏ) - (μ/4))].

The first term (1/ℏ ω) is named Thomas Fermi term and it is apparently the mean of the expression. μ is the Maslov index, which for usual systems is 2 (the number of "regular" turning points in an orbit). The remnant will be called  (Subscript[δg, SC]) below.

As the infinite summation is usually unfeasible one usually consider a "smoothing" of these expressions, suppressing higher harmonics. Let γ the smoothing parameter. The δ functions in  () are replaced by resonance terms, while the semiclassical approximation is smoothed as shown below

g[E] → (∑_(n=1)^∞)((γ/π)/((((E-(Subscript[E, n])))^2)+(γ^2)));
(Subscript[g, SC])[E] → (Subscript[g, TF])(E) + (2/(ℏω)) (∑_(k=1)^∞)(e^(-γ (Subscript[T, k])))Cos[2π k ((I/ℏ) - (μ/4))];     (Subscript[T, k])=((2πk)/ω) .

(Subscript[T, k]) is the period for an orbit covered k times.

We want to check the approximation and the ability of semiclassical approximation to reproduce the spectrum for the systems listed at the beginning of the notebook.

In the following we will need some numerical energy eigenvalues. We will perform a variational computation using the function hermiteMatrices to generate the Hamiltonian terms. The computational details are discussed in notebook [NB-10.10]. The function hermiteFunctions for fixed parity compute the matrices (p^2), (x^2), (x^4) from which the Hamiltonian can be written.

Code for variational computations

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

Harmonic oscillator

For an harmonic oscillator  as it is well known (see chapter on WKB)

E = I ω;   I = (n+(1/2))ℏ

eq.() takes the form

(Subscript[g, SC])[E] = (1/(ℏω))+(2/(ℏω)) (∑_(k=1)^∞)Cos[2π k ((E/(ℏω)) - (1/2))]

In the following we will use ℏ = 1 = ω. We have then, putting iA for the action integral

iA[ene_]:= ene;
gTF[ene_]:= 1;
dgSC[ene_,kmax_]:= 2 Tr[Cos[2π (iA[ene]-(1/2)) Range[kmax]]];

We can simply see a "clustering" around the exact quantum levels (this must be true as the semiclassical approximation is exact for the harmonic oscillator spectrum)

Plot[gTF[x]+dgSC[x,10],{x,0,5},PlotRange-> All]

A quantitative comparison with the exact result which can be generalized to other systems is obtained by the smoothing procedure

The theoretical (exact) energy levels and the smoothing spectral density are

eneTH = 0.5+Range[0,50];
gThSmooth[ene_,eneTH_,γ_]:= Tr[((γ/π)/((((ene-eneTH))^2)+(γ^2)))];

while the semiclassical smoothed density is

gSCSmooth[ene_,kmax_,γ_]:= gTF[ene]+2 Tr[Exp[-γ 2π Range[kmax]] Cos[2π (iA[ene]-(1/2)) Range[kmax]]]

Here is the comparison for γ = 0.1 and kmax = 50. The dashed line represents exact results.

γ = 0.1; kmax = 50;
Plot[{gThSmooth[x,eneTH,γ],gSCSmooth[x,kmax,γ]},{x,0,5},PlotStyle-> {{Red,Dashing[0.015]},Black},PlotLabel-> "γ = "<>ToString[γ],
AxesLabel-> {"E","g(E)"}]

The two density are almost indistinguishable.

The pure quartic oscillator

We consider the potential

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

for g = 1.

g = 1;

Variational results for energy levels

Using a notation familiar from the chapter on variational computations, see as an instance [NB-10.10], we can easily extract the first 50 eigenvalues for the problem:

dim = 200;

α=x/.FindRoot[(2dim+1)(x-(0/(x^3)))-((3*g)/(x^5))(2(dim^2)+2dim+1),{x,(((3g))^(1/3))}];
{p2,x2,x4,x1} = hermiteMatrices[dim];
H0 = (1/2)(α^2)p2 + (0/(2(α^2)))x2 + (g/(2(α^4)))x4;
eneTH= Sort[Eigenvalues[H0]][[Range[50]]];

Semiclassical integrals

The action is

I = (2/(2π))(∫_(-a)^a) dx(Sqrt[2 E - g x^4]) ;    a = (((2E/g))^(1/4)) ;

The result is

I = (((2^(3/4)) (ene^(3/4)) Γ[(5/4)])/((Sqrt[π]) Γ[(7/4)] (g^(1/4))))

as can be easily checked.The frequency is by definition

ω = ((d E)/(d I)) = 1/(dI/dE);

We can write

iA[ene_]:= (((2^(3/4)) (ene^(3/4)) Gamma[(5/4)])/((Sqrt[π]) Gamma[(7/4)] (g^(1/4))));
ω[ene_]:=1/( (3/4)iA[ene]/ene);

Computation of the spectral densities

We take γ = 0.2

γ = 0.2;

The numerical spectral density is always defined by the function gThSmooth, the semiclassical density is given by

gTF[ene_]:= (1/(ω[ene]));
gSCSmooth[ene_,kmax_,γ_]:= gTF[ene]+(2/(ω[ene])) Tr[Exp[-γ(( 2π* Range[kmax])/(ω[ene]))] Cos[2π (iA[ene]-(1/2)) Range[kmax]]]

kmax = 50;
Plot[{gThSmooth[x,eneTH,γ],gSCSmooth[x,kmax,γ]},{x,0.2,20},PlotStyle-> {{Red,Dashing[0.015]},Black},PlotLabel-> "γ = "<>ToString[γ],
AxesLabel-> {"E","g(E)"}]

The agreement is again impressive. This is another way to see that for large quantum numbers the semiclassical approximation works.

The anharmonic oscillator

We take the model

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

with μ = 1 and g = 1:

g=1.; μ = 1.;

Let us proceed as in the previous example

Variational results for energy levels

Using a notation familiar from the chapter on variational computations, see as an instance [NB-10.10], we can easily extract the first 50 eigenvalues for the problem:

dim = 200; μ = 1.;

α=x/.FindRoot[(2dim+1)(x-(μ/(x^3)))-((3*g)/(x^5))(2(dim^2)+2dim+1),{x,(((3g))^(1/3))}];
{p2,x2,x4,x1} = hermiteMatrices[dim];
H0 = (1/2)(α^2)p2 + (μ/(2(α^2)))x2 + (g/(2(α^4)))x4;
eneTH= Sort[Eigenvalues[H0]][[Range[50]]];

Semiclassical integrals

The action and the frequency for this model have been computed in the problems of chapter 10.  Here we quote the result, directly as a Mathematica function

iA[ene_]:= Module[{aL,aR,Jval},

aL = (Sqrt[ ((-1 +(Sqrt[ 1+8 ene  g]))/(2g))]);aR = (Sqrt[ ((1 +(Sqrt[ 1+8 ene  g]))/(2g))]);
Jval = (2/(3π))(Sqrt[g])(Sqrt[(aR^2)+(aL^2)])((aR^2) EllipticK[((aL^2)/((aR^2)+(aL^2)))] - ((aR^2)-(aL^2)) EllipticE[((aL^2)/((aR^2)+(aL^2)))]); Jval]

ω[ene_]:= Module[{aL,aR,ris ,TW},
aR =(Sqrt[ ((1+ Sqrt[1+8ene g])/(2g))]);aL = (Sqrt[ ((-1 + Sqrt[1+8ene g])/(2g))]);
TW = (4/(Sqrt[g]))(1/(Sqrt[(aR^2)+(aL^2)]))EllipticK[((aL^2)/((aL^2)+(aR^2)))];
ris = 2π/TW]

Computation of the spectral densities

We have just to copy the above lines. We take γ = 0.1

γ = 0.1;

The numerical spectral density is always defined by the function gThSmooth, the semiclassical density is given by

gTF[ene_]:= (1/(ω[ene]));
gSCSmooth[ene_,kmax_,γ_]:= gTF[ene]+(2/(ω[ene])) Tr[Exp[-γ(( 2π* Range[kmax])/(ω[ene]))] Cos[2π (iA[ene]-(1/2)) Range[kmax]]]

kmax = 50;
Plot[{gThSmooth[x,eneTH,γ],gSCSmooth[x,kmax,γ]},{x,0.2,20},PlotRange-> All,PlotStyle-> {{Red,Dashing[0.015]},Black},PlotLabel-> "γ = "<>ToString[γ],
AxesLabel-> {"E","g(E)"}]

The agreement is again very good.

The meaning of the Thomas Fermi term

In which sense the Thomas Fermi term is an average distribution is best appreciated by smoothing levels with a relatively big value for γ and comparing the result with the TF term

γ = 1.5;

kmax = 50;
Plot[{gTF[x],gThSmooth[x,eneTH,γ]},{x,0.2,30},PlotRange->{ All,{0,0.5}},PlotStyle-> {{Red,Dashing[0.015]},Black},PlotLabel-> "γ = "<>ToString[γ],
AxesLabel-> {"E","g(E)"}]