Metastable states and  scattering

Quantum Mechanics
                      K. Konishi, G. Paffuti

Study numerically the time evolution of metastable states in the radial potential

V[r]= (Subscript[V, 0]) (r^2)Exp[-r];   (Subscript[V, 0])= 7.5 .

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

Introduction

In this notebook we make a simulation of a simplyfied scattering experiment. We study the time evolution of an s-wave packet impinging onto a central potential

V[r] = (Subscript[V, 0])(r^2)Exp[-r]

We shall see that the incident wave is partly "captured" by the potential forming a metastable state which decays slowly after its formation. A Fourier analysis of the trapped wave reveals the energy of the formed state. The results are compared with analogous results obtained with WKB, Gamow Siegert equation and phase shifts analysis.

A similar simulation has been performed in notebook [NB-3.5-Evolution-tunneling.nb] for the simple case of a square well potential. We put the system in a large box for convenience. Then we approximate time evolution operator for small time steps Δt by

U[Δt] ∼ 1 - i H Δt ∼ ((1 - i H ((Δt)/2))/(1 + i H ((Δt)/2)))

This approximation preserve unitarity. The Hamiltonian matrix H is computed using a finite approximation of Hilbert space spanned by a Fourier basis. In the box of length L a basis is given by

(Subscript[f, n])(x) = (Sqrt[(2/L)])Sin[((n π x)/L)]

Matrix elements of H are easily computed using the Mathematica functions FourierCos and FourierSin. The procedure has been analysed in notebook [*] and the functions developed there are reported below in the section Functions.

Functions

The function getMatV needs as input the values of the potential on the extended grid, i.e. with the extrema, it gives as output the matrix ⟨n|V|k⟩

getMatV[Vx_]:=  Module[{mat,fV,nV,fV1,dim},
dim = Length[Vx]-2;
mat = ConstantArray[0.,{dim,dim}];
nV = Length[Vx]-1;
fV = FourierDCT[Vx,1];
fV1 = Drop[fV,-1];
fV = Join[fV,Reverse[fV1]];
Table[mat[[i,j]] = fV[[i-j+1]] - fV[[i+j+1]],{i,1,dim},{j,1,i-1}];
mat = mat+Transpose[mat];
Table[mat[[i,i]] = fV[[1]] - fV[[2 i+1]],{i,1,dim}];
(1/(Sqrt[2 nV]))mat]

Time evolution

We shall study the time evolution in a potential of the form

V0 = 7.5;
V[r_]:= V0 (r^2)Exp[-r];

We will always use the natural units, ℏ = m = 1.

The discretization is done in a finite box of length L. In this box the first dim components of a Fourier basis (i.e. free particle's states)

(Subscript[f, n])(x) = (Sqrt[(2/L)])Sin[((n π x)/L)]

are used to build an Hamiltonian matrix. The use of Fast Fourier Transform (FFT) in the computation of matrix elements requires to define a grid in coordinate space. As explained in notebook [*] the grid must contain dim internal points, and the function getMatV which computes matrix elements of the potential require the value of V[r] on these points and on the extrema of the grid, here r=0 and r=L.

Definition of the grid

L = 60; dim = 512; (* All internal points *)
Npt = dim+1;
dx = L/Npt;
xptTot = L/(Npt)Range[0,Npt]; (* N+1 points *)
xgrid = xptTot[[2;;dim+1]]; (* the internal points of the grid *)
Vx = V[xptTot];

The Hamiltonian

Once the grid has been assigned the Hamiltonian is easily computed, the potential part is given by getMatV, the kinetic part is (k^2)/2 where k are momenta in the box. k = π n/L, n = 1,...

p2 =N[DiagonalMatrix[((( π/L*Range[1,dim]))^2)]];
matV = getMatV[Vx];
ham =  (1/2)p2 + matV;

Infinitesimal time evolution

We choose a small time step dt and compute the approximate time evolution operator

dt= 0.02;

U1 = IdentityMatrix[dim] - I dt/2* ham;
U2 = Inverse[IdentityMatrix[dim] + I dt/2* ham];
dU = U1.U2;

Incident wave

We know that this system has a metastable state with (Subscript[E, 0])∼3.42. We write a gaussian packet centered in (Subscript[k, 0])= (Sqrt[2 (Subscript[E, 0])])of width Δ. Its initial position, xP, is taken outside the range of the potential. Δ must be chosen with a compromize between energy definition and spreading in r: we are working in a finite box then we will have reflected waves on the unphysical right boundary r = L. Δx must be small enough that these reflections can be neglected during time evolution. We will use a parameter (Δ^2)= 1/3 in k-space, this will correspond in r-space to a gaussian Exp[-3(x^2)].

x0P = 15;
Δ2 = 2./3;
Lrange= 18.;
normψ = 10;

ψ0[ene_,x_]:=normψ* Exp[- I (Sqrt[2. ene]) x - ((Δ2(((x-x0P))^2))/2)]

This is the initial condition. The wave function has been arbitrarily normalized to 10 for graphycal reasons.

Plot[{V[x],Abs[ψ0[(Sqrt[2 3.426]),x]]},{x,0,Lrange},PlotStyle-> {Red,Black},PlotRange-> {All,{0,1.1normψ}}]

Wave function

The passage from k - space to x - space, i.e. the Fourier Transform, is made with FourierDST (Fourier Discrete Sine Transform). The correct normalization to pass from k-space to x - space is (see notebook [*])

ψ[x]=FourierDST[ψk,1]/(Sqrt[dx]);

Where ψk are the Fourier components of ψ. As in this problem the overall normalization is unimportant we omit the prefactor and use simply FourierDST[ψk,1]. This function is its own inverse, i.e. in this normalizartion ψk = FourierDST[ψx,1], where ψx is computed on the internal points of the grid.

Here is how it works :

ψx = ψ0[ 3.426,xgrid];
ψk = FourierDST[ψx,1];

The time evolution

eneResonance = 3.42639;
Γ = 0.0255;

A spread of  δk in k correponds approximatively to a spread of

δE = ∼ (1/2)((((Subscript[k, 0])+ δk))^2)- (1/2)(Subscript[k, 0]^2)∼ (Subscript[k, 0])δk = (Sqrt[2 E Δ2])

in energy.

ene = eneResonance;

numFrames = 100;
dtStepsPerFrame = 10;

getPlot[ene_,numSteps_]:= Module[{},
ψxIn = ψ0[ ene,xgrid];
ψk=FourierDST[ψxIn,1];
 Table[Do[ψk = dU.ψk;,{dtStepsPerFrame}];
ψx= FourierDST[ψk,1];
ListPlot[{Transpose[{xgrid,Abs[ψx]}],Transpose[{xgrid,V[xgrid]}]},Filling-> Axis,
PlotLabel-> j,
Joined-> True,PlotRange-> {{0,L},{0,1.1normψ}}],
{j,1,numSteps}]];

evolutionPlot= getPlot[eneResonance,numFrames];
ListAnimate[evolutionPlot,AnimationRunning->False]

As E shifts from the value at resonance a smaller and smaller part of the incident wave is trapped (proportional to the Fourier component relative at the resonance, see below).

δE = (Sqrt[2 eneResonance Δ2]);

evolutionPlot2= getPlot[Max[eneResonance+2δE,0],Floor[2 numFrames/3]];
ListAnimate[evolutionPlot2,AnimationRunning->False]

Remove[evolutionPlot2]

Spectral analysis

The capture of a resonant state in the scattering can be vividly shown in k-space: if we make the Fourier transform of the reflected wave there must be a "hole" in correspondence of the modes trapped in the potential.

Let us consider again the time evolution for E ∼ (Subscript[E, 0]). We have to choose a time large enough to have a clear distinction between reflected wave and metastable state, but not too long in order to avoid effects due to the reflection on the right boundary. Looking at the first animation one sees that the frame 75 is a good candidate, here we reproduce the relative plot

evolutionPlot[[75]]

We compute again time evolution up to this step

ene = eneResonance;
maxStep = 75;
getPlot[eneResonance,maxStep];

The function getPlot has computed the evolution of ψ and its last value on the grid is ψx, whose absolute value is plotted above. The grid has dim=512 points extending up to x=L, we want to cut points up to x∼8 to select reflected wave. We can drop these first values putting zero for ψ in these places.

ψRflected = PadLeft[Drop[ψx,Floor[6*dim/L]],dim,0.];

Then we make a Fourier transform to extract modes distribution

ψkRflected = FourierDST[ψRflected,1];

kList =N@( π/L Range[dim]);
eneList = (kList^2)/2;
k0 = (Sqrt[2 eneResonance])

spectrum = (L/π)((Abs[ψkRflected])^2)/Tr[((Abs[ψkRflected])^2)];

In the plot below the red dashed line is the theoretical position of the resonance

ListPlot[Transpose[{kList,spectrum}],PlotRange-> {{0,10},{0,1}},Joined-> True,
AxesLabel-> {"k","ψ[k]^2"},
Epilog-> {Dashing[0.02],Red,Line[{{k0,0},{k0,1}}]}]

To make a quantitative analysis the reader can try to work with a larger number of points, dim = 1024 or more, and fit the difference of initial spectrum and reflected wave spectrum by a Lorentzian (translated as a distribution in k).