Spectral density in d=2

Quantum Mechanics
K. Konishi, G. Paffuti

Use the Berry-Tabor formula to obtain the semiclassical approximation for the spectral density for the potential

V[x,y] = ((x^4)/4)+ g ((y^4)/4).

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

Introduction

Summary of the theory

In this notebook we shall use the units,  ℏ = 1 = m.

In the text we reported a semiclassical formula for the spectral density for separable multidimensional system due to Berry e Tabor. In two dimensions is

(Subscript[g, SC])[E] = (Subscript[g, TF])[E]+ ((((2π))^(3/2))/π)(∑_Γ)((Cos[2π M I - σ (π/2)+ 3(π/4)])/(Sqrt[T det[K] (Ω (K^(-1))Ω)])).

The sum runs over the lattice of periodic orbits. (Subscript[Ω, i])are the frequencies and the vector of integers M is defined by

Ω = ((2π)/T) M = ((2π)/T)((Subscript[M, 1]),(Subscript[M, 2])).

The rational nature of (Subscript[M, 1])/(Subscript[M, 2])ensures periodicity.

K is the Hessian of the Hamiltonian matrix, (Subscript[K, ij])= (∂^2)H/∂(Subscript[J, i])∂(Subscript[J, j]), (Subscript[J, i])being the classical action integral. Ω are the corresponding frequencies and T the period of the orbits. (Subscript[g, TF])is the semiclassical Thomas Fermi density, which in two dimension is

(Subscript[g, TF])(E) = (1/(2π))∫Θ(E - V(r))dr.

The limitation of the sum in () to periodic orbits is non trivial and reflect the importance of these orbits for the quantization problem. We note en passant that quantum number selected by WKB conditions (Subscript[J, i])= ((Subscript[n, i])+ 1/2) does not select in general periodic orbits.

We shall apply () to a simple potential

V[x,y] = ((x^4)/4)+ g ((y^4)/4).

The classical actions

The system is separable, the "energies" (Subscript[E, x]) and (Subscript[E, y]) relative to x and y coordinates are separately conserved. The two action variable in x and y coordinate are (see [NB-13.7] )

(Subscript[J, 1])=(Subscript[k, 1]) (Subscript[E, x]) ;              (Subscript[J, 2])=(Subscript[k, 1]) (Subscript[E, y]) (g^(-1/4));           (Subscript[k, 1]) =(((2^(3/4)) Gamma[(5/4)])/((Sqrt[π]) Gamma[(7/4)]));

The energy of the system is E = (Subscript[E, x])+ (Subscript[E, y]) or inverting the relations ():

E = k ( ((Subscript[J, 1])^(4/3))+ (g^(1/3))((Subscript[J, 2])^(4/3)));            k = ((((((Sqrt[π]) Gamma[(7/4)])/((2^(3/4))  Gamma[(5/4)]))))^(4/3)).

The frequencies are easily computed, and the periodicity condition is

Ω = (4/3)k{((Subscript[J, 1])^(1/3)),(g^(1/3))((Subscript[J, 2])^(1/3))} = ((2π)/T){(Subscript[M, 1]),(Subscript[M, 2])}.

We can express action integrals through M

{(Subscript[J, 1]),(Subscript[J, 2])}= ((((3/(4k))((2π)/T)))^3){((Subscript[M, 1])^3),(1/g)((Subscript[M, 2])^3)};

and inserting this expression in the definition of E[J]

E = k ((((3/(4k))((2π)/T)))^4)(((Subscript[M, 1])^4)+((g^(1/3))/(g^(4/3)))((Subscript[M, 2])^4))= k ((((3/(4k))((2π)/T)))^4)(((Subscript[M, 1])^4)+(1/g)((Subscript[M, 2])^4));
T = (k^(1/4))((3/(4k))2π)(((((Subscript[M, 1])^4)+(1/g)((Subscript[M, 2])^4)))^(1/4))/(E^(1/4)).

It is easy to show that for the matrix K, denoted below with Kmat

Kmat = (4/9)k(((Subscript[J, 1])^(-2/3)),(g^(1/3))((Subscript[J, 2])^(-2/3)));
Ω (Kmat^(-1))Ω = ((((4/3)k))^2)(9/(4k))(((Subscript[J, 1])^(4/3))+(g^(1/3))((Subscript[J, 2])^(4/3))) = 4 k (((Subscript[J, 1])^(4/3))+(g^(1/3))((Subscript[J, 2])^(4/3))) = 4 E .

We see that given a frequency lattice we can compute everything for a fixed E, the variable quantity in (Subscript[g, SC])[E].

Comparison with numerical results

To compare a finite truncation of () to numerical data we use a smoothing procedure, see notebook [NB-13.7]

(Subscript[g, SC])[E] = (Subscript[g, TF])[E]+ ((((2π))^(3/2))/π)(Underscript[∑ , Γ])((Cos[2π M I - σ (π/2)+ 3(π/4)])/(Sqrt[T det[K] (Ω (K^(-1))Ω)])) (e^(-γ T))

This is compared with the corresponding smoothed exact spectral density

g[E] → (∑_(n=1)^∞)((γ/π)/((((E-(Subscript[E, n])))^2)+(γ^2)));

In the following we will use the eigenvalues computed with a variational procedure and use the usual function hermiteMatrices[] to write Hamiltonian matrix elements, see notebooks  [NB-10.10], [NB-13.7] .

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

The model

We use

g=1.5; γ =0.2;

for the coupling constant and the smoothing parameter respectively.

The variational results

The first part of the computation is standard. The system is considered as a sum of a subsystem with coupling g and a subsystem with coupling 1. We first compute the eigenvalues for the subsystem with coupling g=1.5

μ  = 0 ;    
dim = 100;
{p2,x2,x4,x1} = hermiteMatrices[dim];

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

then the eigenvalues for the subsystem with g=1 :

g1=1.;
α=x/.FindRoot[(2dim+1)(x-(μ/(x^3)))-((3g1)/(x^5))(2(dim^2)+2dim+1),{x,(((3g1))^(1/3))}];
H0 = (1/2)(α^2)p2 + (μ/(2(α^2)))x2 + (g1/(2(α^4)))x4;
e2= Sort[Eigenvalues[H0]];

The function Outer[Plus,e1,e2] sums in all possible ways the pair of eigenvalues computing the spectrum of the composite system. See Mathematica help for Outer.

eneTH=Sort[Flatten[Outer[Plus,e1,e2]]];

The smoothed density is given by

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

The semiclassical computation

Thomas Fermi density

The definition can be written

(Subscript[g, TF])(E) = (1/(2π))∫Θ(E - V(r))dr = (1/(2π))(∫_(-a)^a)dx (∫_(-b)^b)dy ;
a = (((2 E))^(1/4));  b = (((((2 E - (x^4))/g)))^(1/4));

The result is

(1/(2π))Integrate[1,{x,- (((2 ene))^(1/4)),(((2 ene))^(1/4))},
{y,-(((((2 ene - x^4)/gC)))^(1/4)),(((((2 ene - x^4)/gC)))^(1/4))},
Assumptions-> {ene > 0,gC>0}]

gTF[ene_]:=(((Sqrt[ene]) ((Gamma[(1/4)])^2))/(2 (Sqrt[2]) (g^(1/4)) (π^(3/2))))

Periodic orbits

We first fix a maximum value for the lattice integers

orbMax = 10;

Then we generate the lattice using again the function Outer. With argument List it generates all pairs  ((Subscript[M, 1]),(Subscript[M, 2]))

mV =Flatten[ Outer[List,Range[orbMax],Range[orbMax]],1];

Then we define a function which compute the oscillating (smoothed) part of the density

k = N[((((((Sqrt[π]) Gamma[(7/4)])/((2^(3/4))  Gamma[(5/4)]))))^(4/3))];

dgSC[ene_]:= Module[{T,Ω,J,detK,Km1,k2,ris,ris1,phase,σ},
(* period *)
T=(k^(1/4))((3/(4k))2π)(ene^(-1/4))(((Tr[{1,(1/g)}*#1^4]&/@mV))^(1/4));
(* frequencies: Ω = 2 π/T (M1,M2) *)
Ω=((2π)/T)mV;
(* action integrals *)
J = ((((3/(4k))((2π)/T)))^3)*({1,(1/g)}*#1&/@ (mV^3));
(* computation of (Sqrt[T det[K] (Ω (K^(-1))Ω)]) *)
σ = 2+({2,2}.#1&/@mV);
detK=Tr[#1,Times]&/@((4/9)k*{1,(g^(1/3))}*#1&/@(J^(-2/3)));
k2 = 4 ene; 
phase = (2π*Tr[#1]&/@ (mV*J)) -(π/2)σ + ((3π)/4);
ris1 = Cos[phase]/(Sqrt[T*detK*k2]);
ris1 = Exp[-γ T]*ris1;
ris = (1/π)(((2π))^(3/2))*Tr[ris1]]

Comparison

Here we plot the two densities

Plot[{gThSmooth[x,eneTH[[Range[50]]],γ],gTF[x] + dgSC[x]},{x,0,20},PlotStyle-> {{Red,Dashing[0.015]},Black},PlotLabel-> "γ = "<>ToString[γ],
AxesLabel-> {"E","g(E)"}]

The details depend on the smoothing parameter but it is evident that the truncated sum over periodic orbits gives quite a good approximation to the quantum spectrum.