Decay

Quantum Mechanics
 K. Konishi, G. Paffuti

Study numerically the evolution problem for the potential

V[x] = ({{∞;,      x< 0}, {0;,       0≤x≤a}, {(Subscript[V, 0]);,      x>a}}).

At t>0 the potential V[x] is taken to vanish for x > 2a. Study the evolution of the metastable state and its decay law.

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

The model

The model has been studied in problem [prob.11]. Consider a particle in a potential well of height (Subscript[V, 0]). Region x < 0 is forbidden by an infinite potential barrier. A particle is confined in the well, in a stationary state. At t=0 the barrier after x = a+δa disappears and the particle can escape via tunnel effect. The situation is depicted below.

({{, }})

We wish to study the time evolution of the "metastable state" and compare the rate of the decay with the theoretical result.

Bound states

Bound states are found by imposing matching conditions at x=a on the wave function

ψ[x] ={ ({{A Sin[k x],    , 0≤ x ≤ a}, {B Exp[-β x],  , x> a}})   ;        with      k = (Sqrt[((2m E)/(ℏ^2))]);  β = (Sqrt[((2m ((Subscript[V, 0])-E))/(ℏ^2))]) .

The problem easily is easily reduced to the computation of the intersections of the two curves

η = - ξ Cot[ξ];    (ξ^2)+ (η^2) = ((2m (Subscript[V, 0])(a^2))/(ℏ^2));             with   ξ=ka;η=βa .

Plot[-ξ Cot[ξ],{ξ,0.001,6},AspectRatio-> Automatic,PlotRange-> {{0,6},{0,5}},
AxesLabel-> {"ξ","η"},Epilog-> Circle[{0,0},2]]

It is clear that we will have at least a solution only for

(π/2) < (Sqrt[((2m (Subscript[V, 0])(a^2))/(ℏ^2))])

We shall choose parameters in such a way that only one bound state will be present (the above parameters must be less than 3π/2 )

N.B.  The units m = ℏ = a = 1 will be chosen below.

a=1.;

Vmin= (((π/2))^2)/(2 (a^2));

a=1; da =  0.6a;
V0 = 8Vmin;

The potential is defined as (we will set the attibute Listable to the function to allow a simple way to work with vectors):

V[x_] := If[x < a || x > a + da, 0, V0];
SetAttributes[V, Listable];

You can use Mathematica  function Piecewise, or an unstructon Which[] instead of the If[].

This is the equation for determining k for the bound state :

kState=Re[ ξ/.FindRoot[-ξ Cot[ξ]== (Sqrt[(2V0)-(ξ^2)]),{ξ,π/2,π}]]

Here we construct the initial wave function

eneState= (kState^2)/2;
β = (Sqrt[2(V0-eneState)]);

ψIn[x_]:= If[x≤ a,Sin[kState*x],Sin[kState*a]Exp[-β( x-a)]];
SetAttributes[ψIn,Listable];

normψIn = 1/Sqrt[NIntegrate[((Abs[ψIn[x]])^2),{x,0,∞}]]

A pictorial view is presented in the graph (ψ has been rescaled for graphical reasons)

barrier[{a1_,b1_},{a2_,b2_}]:= With[{d=200,dL= b1-a1},
Table[{ColorData["GrayTones"][n/d/1.1],
Rectangle[{a1+dL*n/d,a2},{a1+dL*(n+1)/d,b2}]},{n,0,d}]];

Plot[5*ψIn[x]*normψIn,{x,0,5a},Axes-> None,PlotRange-> {{-0.3a,5a},{-0.1,1.1V0}},
Epilog-> {Line[{{0,0},{a,0},{a,V0},{5a,V0}}],
barrier[{0,-1/6},{0,2V0}]}]

Time evolution

The system is expected to decay at a rate (see problems ) that can be estimated in

W = (dP/dt) = (v/(2a))D ;    D =  ((4 (k^2) (β^2))/(4 (k^2) (β^2)+((((k^2)+(β^2)))^2) ((Sinh[β δa])^2))).

To study time evolution we adopt the simple method of discretization,  i.e. we use a spatial grid of spacing h and approximate infinitesimal time evolution operator by

dU ≃ ((1 - i (1/2) H Δt)/(1 + i (1/2) H Δt)).

The function getdU has been already used in problem [], the only modification here is the we define the starting point of the grid at x=0. The output will consist of grid points and infinitesiamal evolution operator, the input require the infinitesimal time step, the potential function, the number of points of the grid and a length L, which is the length of a big box in which we imagine the systemis plonged. In the solution we must take care of choosing times in such a way to avoid reflections on the x = L boundary.

getdU[Δt_,V_,nP_,L_]:= Module[
{D2,indDiag,matV,h,x,Ham,dU,U2,U1,vecV},
(* definition of the grid *)
h = N[L/(nP+1)]; 
x =  h* Range[nP];
indDiag = Table[{i,i},{i,1,nP}];
(* Kinetic term *)
D2 = SparseArray[{{i_,i_} ->  -2,{i_,j_}/;j-i==1-> 1,
{i_,j_}/; i-j== 1->  1},{nP,nP}];
(* potential term and Hamiltonian *)
(* vecV = V[#1]&/@ x; *)
matV = SparseArray[indDiag-> V[x]];
Ham = -1/(2 (h^2)) D2 + matV;
(* evolution operator *)
U2 = IdentityMatrix[nP] - I Δt/2 Ham;
(* U1 = IdentityMatrix[nP]; *)
U1 =Inverse[ IdentityMatrix[nP] +  I Δt/2 Ham]; 
dU = U2.U1;
{x,dU}]

Here we choose the grid. We take L = 40 and a time interval Δt = 0.025. We produce values for ψ for a certain number of "frames", the interval between two frames has been chosen (dtStepsPerFrame × Δt). The partial results are stored in a matrix matPsi whose rows contains the value of ψ at the grid points.

The reader can change these parameters, maybe worsening the accuracy: the procedure is a bit time-consuming with the proposed values.

L = 40; nP = 1000; 
Δt = 0.025; dtStepsPerFrame = 4; tmax=15;
numFrames = Floor[tmax/(dtStepsPerFrame*Δt)];

The following steps have been repeatedly used  in others notebooks.

{xpt,dU} = getdU[Δt,V,nP,L];
ψ = normψIn*ψIn[xpt];
maxψ = Max[Abs[ψ]];
matPsi = Table[Do[ψ = dU.ψ;,{dtStepsPerFrame}];ψ,{numFrames}];

ListAnimate[ Table[
ListPlot[Transpose[{xpt,Abs[matPsi[[j]]]^2}],Filling-> Axis,Joined-> True,PlotRange-> {{0,L},{0,1.1((maxψ)^2)}}],
{j,1,numFrames}],AnimationRunning->False]

The decay of a metastable state

We will use the precedent computation to study the probability of finding the initial state in the evolving ψ,

P[t] = ((|⟨(Subscript[ψ, 0])|ψ[t]⟩|)^2)

An exponential decay correspond to

P[t] = Exp[-Γ t];

τ = 1/Γ  is the mean life of the state.

The semiclassical approximation given above would predict :

gammaTheo[k_,V0_,δa_]:= Module[{ene,β,dF},
ene = k^2/2; β = (Sqrt[2(V0-ene)]);
dF =((4 (k^2) (β^2))/(4 (k^2) (β^2)+((((k^2)+(β^2)))^2) ((Sinh[β δa])^2)));(k/(2a))dF]

γTheo= gammaTheo[kState,V0,da]

Within our approximations it is faster and it is accurate  enough to compute the integrals needed for P[t] using trapezoidal rule:

(∫_a^b)f[x]𝕕x ≃ h( (1/2)f[(Subscript[x, 0])] + f[(Subscript[x, 1])] + …  + f[(Subscript[x, N-1])]+ (1/2)f[(Subscript[x, N])])

Our functions at end points are (approximately) zero, so the previous formula reduces to the sum of values times the space grid step h. As an example the norm of initial wave function is, in this approximation:

prob0 = (xpt[[2]]-xpt[[1]])Total[((ψIn[xpt])^2)];

To be compared with the coputed value 1/((normψIn)^2)

{prob0,1/((normψIn)^2)}

The list of probabilities as a function of time is then  (we prepend the normalized value at t=0, i.e. 1)

prob = Prepend[Table[((((xpt[[2]]-xpt[[1]])*Total[Abs[matPsi[[j,1;;All]]*normψIn*ψIn[xpt]]]))^2),{j,1,Length[matPsi]}],1.];

Let us fit our results with an exponential

times = Range[0,numFrames]Δt*dtStepsPerFrame;

data=Transpose[{times,prob}];
model = a1*Exp[-γ t];
fit=FindFit[data,model,{a1,γ},t]

modelf=Function[{t},Evaluate[model/.fit]];
Plot[modelf[t],{t,0,tmax},AxesLabel-> {"t","P[t]"},
Epilog->Map[Point,data[[1;; ;;5]]]]

Concluding remarks and exercises

Small time behaviour

There is a very good exponential behavior, even if at small times there is a small discrepancy:

Plot[modelf[t],{t,0,2},AxesLabel-> {"t","P[t]"},
Epilog->Map[Point,data[[1;;20]]]]

We will see in the study of perturbation theory that in effect for very small times there must be a quadratic dipendence on time.

See the related remarks on Quantum Zeno effect in the main Text, in Section 13.2.2.

Result on γ

The agreement with semiclassical formula (.) is qualitatively correct, but as we will see with some details studying WKB approximation, semiclassical results work better for higher excited states.

The result do not depend on our numerical approximation for integrals, the reader can easily construct a table for probabilities using intepolating functions, nothing qualitatively new happens.

Later in this book we will see that γ can be computed by solving a time independent Schrödinger equation, with the parameters used in this example, (Subscript[V, 0])= 8 (π^2)/(8 (a^2)), the result wil be γ = 0.043.

Observables

The reader can try some other observables. For instance instead of the projection on the initial state a more reasonable observable, in this case, is the probability of finding the particle inside the region [0, a + δa].

As an exercise the reader can follow the evolution for a "non metastable" state, i.e. a state which do not correspond to a bound state at t=0. As an example for k0 = π/a one find for the surviving possibility

Can the reader explain this behavior? Why does an exponential behavior seem to reappear at long times? What is the behavior of the probability of the particle being in [0, a + δa] ?

NDSolve

It is worth to check the results with a short time computation using NDSolve. The calculation take some time, it is convenient to work with smaller L and tmax, checking on the animation that reflection on x=L are negligible. One possible chuck of code is (L=10, tmax = 3):

equation[L_]:= {i(∂_t)u[x,t]==-1/2(∂_(x,x))u[x,t]+V[x]u[x,t] ,
u[x,0]==ψIn[x],u[L,t]==0, u[0,t]== 0}

ψ = u /. First[NDSolve[equation[10], u, {x, 0, 10}, {t, 0., 3},
      PrecisionGoal -> 2, AccuracyGoal -> 2 ]]; // Timing

The reader can explore different options in NDSolve, see Mathematica help for the various possibilities.