Density of states and scattering phases 2

Quantum Mechanics
 K. Konishi, G. Paffuti

Compute the resonance parameters by scattering phases in the one dimensional potentials (x^2)/2(1+(gx^6)) and (x^2)/2(1+(gx^4)).  
To compute the scattering phases use both the transmission and reflection coefficients and the regular solution of the Schrödinger equation.

ClearAll["Global`*"];
Off[General::spell, General::spell1];
Off[InterpolatingFunction::"dmval"];
Off[Eigenvalues::"arh"];

Introduction

Summary of the theory

In the text it has been pointed out that the spectral density (at fixed parity in one dimensional problems) is related to the scattering phases by

g[E] = (Subscript[g, 0])[E] + (1/π)((dδ)/dE).

Where (Subscript[g, 0]) is the free particle contribution. In particular near a resonance

g[E]≃ (1/π) ((Γ/2)/((Γ^2)/4 + (((E-E0))^2))).

In this notebook we will compute g[E], i.e. scattering phases, by two different methods.

Transmission and reflection coefficients

In one dimension it is possible to solve directly the scattering problem and extract form the solution the phases. In the text it has been shown that the general form of transmission and reflection coefficients in a parity preserving problem is

(Subscript[A, T])= Cos[θ]Exp[i α];    (Subscript[A, R])= i Sin[θ]Exp[i α].

The scattering phases in even and odd channel are given by

(Subscript[δ, e])= ((α + θ)/2) ;   (Subscript[δ, o]) = ((α - θ)/2).

In this framework we can use the same procedure used in the notebooks [***]. We write a wave function in the right asymptotic zone as a positive exponential and match the solution on the left asymptotic zone with a superposition of incident and reflected waves

left: a Exp[i k x] + b Exp[-i k x];            right:   Exp[i k x]

The coefficients (Subscript[A, T]) and (Subscript[A, R]) follow

(Subscript[A, T])= 1/a ;      (Subscript[A, R])= b/a .

Regular functions

In one dimension the regular functions have an asymptotic behavior Cos[kx + δ] and  Sin[kx + δ] (even / odd functions). This is a good exercise to work for the reader, and we refer to the problems in this chapter for the proof.

In this scheme the computation of the phases is almost identical to the 3-dimensional case. Let  ψEven,  ψOdd the regular solutions with fixed parity. For a smooth potential they satisfy the boundary conditions

ψEven[0]=1;ψEven'[0]=0;               ψOdd[0]=0;ψOdd'[0]=1;

Let us now call fB1 and fB2 two elementary solutions which behave as Cos[k x] and Sin[k x] asymptotically (fB1 and fB2 can be Cos and Sin itself ). Asymptotically we have

ψEven ∼ A Cos[δ] fB1 - A Sin[δ] fB2 ;     ψOdd ∼ A Sin[δ] fB1 + A Cos[δ] fB2 ;

Given two asymptotic values for ψ and computing there fB1 and fB2, the phase δ can be extracted, as in the radial equation for the 3-dimensional case.

Asymptotic solutions

Usually fB1, fB2 are Cos[], and Sin[], but if the potential goes to zero very slowly, let us say as 1/(x^2), it can be convenient to change the basis.

For a potential decreasing as α/2(x^2) the asymptotic equation takes the form

ψ'' + k^2 ψ = α/x^2   ψ;   i.e.:  ψ'' + k^2 ψ - (((ν^2) - 1/4)/(x^2))ψ = 0;     ν = (Sqrt[α+(1/4)])

The general solution is

A (Sqrt[x ])BesselJ[ν,kx]+ B (Sqrt[x ])BesselY[ν,kx]

From the asymptotic  expansions

(Subscript[J, ν])(x)∼(Sqrt[(2/(π x))])Cos[x-ν(π/2)-(π/4)];    (Subscript[Y, ν])(x)∼(Sqrt[(2/(π x))])Sin[x-ν(π/2)-(π/4)];

we see that it is convenient to take as basis function

fB1 = (Sqrt[x ])BesselJ[ν,kx];    fB2 = (Sqrt[x ])BesselY[ν,kx]

These give,  through eq.(), the phase δ, apart for a constant -π/2ν - π/4 which is inessential for δ'[E].

Some utility functions

To smooth the discontinuity introduced by of ArcTan we use the function smoothPhase, introduced in notebooks [NB-11.7], [NB-13.5]:

δ = smoothPhase[phase];

The system () is solved by computing the solution at two values of x and extracting the transmission and reflection coefficients

{AT, AR} = getAtAr[val1, val2, x1, x2, k];

If the regular solution approach is followed we have to solve eq.(). The strategy is always to compute the exact solution in two points and extract δ. We give two almost identical functions, one solve the system for a base Cos, Sin, the other for a base J, Y. In the second case the parameter ⋁ must be supplied

phase = phaseCS[val1, val2, x1, x2, k, parity];

phase = phaseJY[ν,val1, val2, x1, x2, k, parity];

Code of the functions

smoothPhase[phase_]:= Module[{disc,jump,outphase},
disc = 0*phase;
(* This detect the jumps *)
Do[ 
jump = phase[[i]] - phase[[i-1]];
disc[[i]]=Which[jump < -1.7*Pi,1,jump > 1.7*Pi,-1,True,0];,
{i,2,Length[phase]}];
(* This add 2π for each jump *)
outphase = Drop[FoldList[Plus,0,disc],1]*2Pi+phase]

getAtAr[val1_,val2_,x1_,x2_,q_]:= Module[{c11,c12,c21,c22,a1,a2,a,b},
c11 = Exp[I q x1];
c21 = Exp[-I q x1];
c12 = Exp[I q x2];
c22 =  Exp[-I q x2];
{a,b}={a1,a2}/.First[NSolve[{a1*c11 + a2 *c21 == val1,a1*c12 + a2 *c22 == val2},{a1,a2}]]; 
{1./a,b/a}]

phaseCS[val1_,val2_,x1_,x2_,q_,parity_]:= Module[{c11,c12,c21,c22,a1,a2,d1,d2,ph},
c11 = Cos[q x1];
c21 = Sin[q x1];
c12 = Cos[q x2];
c22 = Sin[q x2];
{d1,d2}={a1,a2}/.First[NSolve[{a1*c11 + a2 *c21 == val1,a1*c12 + a2 *c22 == val2},{a1,a2}]];
(*If[parity== 0,ArcTan[d1,-d2],ArcTan[d2,d1]]*)
If[parity==0,ph = Arg[d1 - I d2];,ph = Arg[d2 + I d1];];
ph]

phaseJY[ν_,val1_,val2_,x1_,x2_,q_,parity_]:= Module[{c11,c12,c21,c22,a1,a2,d1,d2,ph},
c11 = (Sqrt[x1]) BesselJ[ν,q x1];
c21 = (Sqrt[x1]) BesselY[ν,q x1];
c12 = (Sqrt[x2]) BesselJ[ν,q x2];
c22 = (Sqrt[x2]) BesselY[ν,q x2];
{d1,d2}={a1,a2}/.First[NSolve[{a1*c11 + a2 *c21 == val1,a1*c12 + a2 *c22 == val2},{a1,a2}]];
(*If[parity== 0,ArcTan[d1,-d2],ArcTan[d2,d1]]*)
If[parity==0,ph = Arg[d1 - I d2];,ph = Arg[d2 + I d1];];
ph]

Example 1

g = 0.02;  
V[x_]:= (1/2)((x^2)/(1+g (x^6)))
 Plot[V[x],{x,-5,5},AxesLabel-> {"x","V[x]"}]

Phases through (Subscript[A, T]) and (Subscript[A, R])

We first choose a far enough point where to compute solutions.

Bx =   100; Δx = 0.1;

Here we define a list of energies, the grid is not uniform.

Emin = 0.001;Emax = 5; 
dE = 0.005; dE2=0.02; 
e1 = Range[Emin,1,dE];
e2 = Range[Max[e1]+dE2,2,dE2];
e3 = Range[Max[e2]+dE2,Emax,2*dE2];
ene = Join[e1,e2,e3];
k = (Sqrt[2 ene]);

We solve the equation and extract the transmission and reflection coefficients with the function getAtAr[]

AT = Table[0.,{Length[k]}]; AR = AT;
Do[
solution=NDSolve[{y''[x] + k[[i]]^2y[x] -2 V[x]y[x] ==0,
y[Bx]== Exp[I k[[i]] Bx],y'[Bx]==I k[[i]]Exp[I k[[i]] Bx] },y,{x,Bx,-Bx}];
val1 = y[-Bx]/.solution;
 val2=y[-Bx+Δx]/.solution;
{AT[[i]],AR[[i]]}= getAtAr[val1,val2,-Bx,-Bx+Δx,k[[i]]];
,
{i,1,Length[ene]}];

The phases are computed with eq.(4).

δEven = smoothPhase[Arg[AT+AR]]/2;
δOdd = smoothPhase[Arg[AT-AR]]/2;

We interpolate to obtain functions δ[E] and define the relative spectral densities by eq.(2)

fδEven = Interpolation[Transpose[{ene,δEven}]];
fδOdd = Interpolation[Transpose[{ene,δOdd}]];
dgEven[x_]:=(1/π)Derivative[1][fδEven][x];
dgOdd[x_]:=(1/π)Derivative[1][fδOdd][x];

GraphicsRow[{Plot[dgEven[x],{x,0.1,0.6},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Even",
AxesLabel-> {"E","δg"}],
Plot[dgOdd[x],{x,0.1,1.5},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Odd",
AxesLabel-> {"E","δg"}]}]

A fit allow an estimation of the even parity resonance:

Clear[a,b,c,d];
ptEne = Range[0.35,0.55,0.005];
ptDelta = dgEven[ptEne];
points = Transpose[{ptEne,ptDelta-ptDelta[[1]]}];
{E0,Γ}/.FindFit[points,b/2/(Γ^2/4+ (x-E0)^2),{b,E0,Γ},x]

Phases through regular solutions

To the same system we can apply the method of regular solutions

Bx=100;

We solve the equations and extract phases by phaseCS function

δEven = Table[0,{i,1,Length[ene]}];
δOdd=δEven;
Do[
solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 1,y'[0]==0 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δEven[[j]]=phaseCS[val1,val2,Bx,Bx-Δx,k[[j]],0];
(* Odd solution *)
 solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 0,y'[0]==1 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δOdd[[j]]=phaseCS[val1,val2,Bx,Bx-Δx,k[[j]],1]; 
,
{j,1,Length[ene]}];

Interpolation, graphs and fit are repeated

fδEven = Interpolation[Transpose[{ene,smoothPhase[δEven]}]];
fδOdd = Interpolation[Transpose[{ene,smoothPhase[δOdd]}]];
dgEven[x_]:=(1/π)Derivative[1][fδEven][x];
dgOdd[x_]:=(1/π)Derivative[1][fδOdd][x];

GraphicsRow[{Plot[dgEven[x],{x,0.1,0.6},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Even",
AxesLabel-> {"E","δg"}],
Plot[dgOdd[x],{x,0.1,1.5},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Odd",
AxesLabel-> {"E","δg"}]}]

Clear[a,b,c,d];
ptEne = Range[0.35,0.55,0.005];
ptDelta = dgEven[ptEne];
points = Transpose[{ptEne,ptDelta-ptDelta[[1]]}];
{E0,Γ}/.FindFit[points,b/2/(Γ^2/4+ (x-E0)^2),{b,E0,Γ},x]

The results are equivalent to the previous ones.

Example 2

As a second example we consider

g = 0.1;xmax = (g^(-1/4));  
V[x_]:= (1/2)((x^2)/(1+g (x^4)))
 Plot[V[x],{x,-5xmax,5xmax},AxesLabel-> {"x","V[x]"}]

Using notebook [NB-13.1] it is easy to compute parameter resonances for g=0.05 and g-0.1:

{Eth05,Γth05}={0.4546,2*0.001449};
{Eth1,Γth1}={0.41507,2*0.0077};

First we apply the same method of regular solutions. We just copy the above instructions,  so no comments are required.

Bx =  50.;  Δx=0.4;

Emin = 0.002;Emax = 3.5; 
dE = 0.001; dE2=0.02; dE3 = 1;
e1 = Range[Emin,0.6,dE];
e2 = Range[Max[e1]+dE2,2,dE2];
e3 = Range[Max[e2]+dE2,Emax,2*dE2];
ene = Join[e1,e2,e3];
k = (Sqrt[2.ene]);

δEven = Table[0,{i,1,Length[ene]}];
δOdd=δEven;
Do[
solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 1,y'[0]==0 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δEven[[j]]=phaseCS[val1,val2,Bx,Bx-Δx,k[[j]],0];
(* Odd solution *)
 solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 0,y'[0]==1 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δOdd[[j]]=phaseCS[val1,val2,Bx,Bx-Δx,k[[j]],1]; 
,
{j,1,Length[ene]}];

fδEven = Interpolation[Transpose[{ene,smoothPhase[δEven]}]];
fδOdd = Interpolation[Transpose[{ene,smoothPhase[δOdd]}]];
dgEven[x_]:=(1/π)Derivative[1][fδEven][x];
dgOdd[x_]:=(1/π)Derivative[1][fδOdd][x];

GraphicsRow[{Plot[dgEven[x],{x,0.1,0.6},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Even",
AxesLabel-> {"E","δg"}],
Plot[dgOdd[x],{x,0.1,1.5},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Odd",
AxesLabel-> {"E","δg"}]}]

We see that the computation is a bit erratic at low energies, where large distances are involved. The situation is worse for lower g.

Bessel

To give a partial stabilization to the above results we repeat the computation using the more refined large distance solution () using Bessel function. The only difference in the code is the use of phaseJY[]. From asymptotic behavior of the potential we get the index of Bessel functions, see eq.():

ν = (Sqrt[(1/g)+(1/4)]);

The rest of the solution is by now standard

δEven = Table[0,{i,1,Length[ene]}];
δOdd=δEven;
Do[
solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 1,y'[0]==0 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δEven[[j]]=phaseJY[ν,val1,val2,Bx,Bx-Δx,k[[j]],0];
(* Odd solution *)
 solution=NDSolve[{y''[x] + k[[j]]^2y[x] -2 V[x]y[x] ==0,
y[0]== 0,y'[0]==1 },y,{x,0,Bx}];
val1 = y[Bx]/.solution;
val2 = y[Bx-Δx]/.solution;
δOdd[[j]]=phaseJY[ν,val1,val2,Bx,Bx-Δx,k[[j]],1]; 
,
{j,1,Length[ene]}];

fδEven = Interpolation[Transpose[{ene,smoothPhase[δEven]}]];
fδOdd = Interpolation[Transpose[{ene,smoothPhase[δOdd]}]];
dgEven[x_]:=(1/π)Derivative[1][fδEven][x];
dgOdd[x_]:=(1/π)Derivative[1][fδOdd][x];

GraphicsRow[{Plot[dgEven[x],{x,0.2,0.6},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Even",
AxesLabel-> {"E","δg"}],
Plot[dgOdd[x],{x,0.1,1.5},PlotRange-> All,ImageSize-> 4 72,PlotLabel-> "Odd",
AxesLabel-> {"E","δg"}]}]

The even resonance can be fitted

Clear[a,b,c,d];
ptEne = Range[0.3,0.6,0.0005];
ptDelta = dgEven[ptEne];
points = Transpose[{ptEne,ptDelta-ptDelta[[1]]}];
{E0,Abs[Γ]}/.FindFit[points,b/2/(Γ^2/4+ (x-E0)^2),{b,E0,Γ},x]

The comparison with results in [NB-13.1] is satisfactory: clearly in the method we are using the subtraction of the background can introduce some error.

{{Eth05,Γth05},{Eth1,Γth1}}