Time evolution

Quantum Mechanics
 K. Konishi, G. Paffuti

Study numerically the transmission of a wave packet through a rectangular barrier. Use both a discretization method and NDSolve to solve the Schrödinger equation.
 Use a Fourier analysis of the solution to get the momentum content of the transmitted wave and compare the result with the theoretical expectations.

Utilities

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

Introduction

In notebook [NB-3.4] we have presented a couple of techniques for the numerical study of time-dependent Schrödinger equation:

i ℏ(∂/(∂t))ψ[t,x] = H ψ[t,x];                 H = -(( (ℏ^2))/(2 m)) ((∂^2)/(∂(x^2))) + V[x];

Two different boundary conditions were described: zero b.c. (equivalent to put the system in an impenetrable box) or periodic boundary conditions. If the box has length L the time evolution can be followed up to times τ where boundary effects come in, example reflection on the boundary in the first case.

In this notebook the same methods are applied to a scattering of a wave packet on a barrier, to give a dynamical picture of tunnel effect. The method used below can be easily adapted to any potential, for comparison with analytical known results we present the case of a rectangular barrier.

In what follows we always use units ℏ = m = 1.

Definition of the system

Potential

The potential barrier is taken to be of rectangular form, with the height (Subscript[V, 0]) placed between aL and aR. It is convenient to express V as an effective momentum, V = (ϵ^2)/ 2. Momentums k > ϵ are above the barrier.

aL = -1.1;aR = 0.9;
ϵ0 = 2.0;  V0 = ϵ0^2/2.; 
V[x_]:=V0*UnitStep[x-aL]UnitStep[aR-x];

Wave packet

We take as incident wave packet a gaussian wave centered in x0 with central momentum k0 and a width parameter dk0

k0 = 0.9 ϵ0;dk0 =  0.15*k0;
x0P = -25;
ψ0[x_]:= (E^(I k0 x- (dk0^2)(((x-x0P))^2)));

Plot[((Abs[ψ0[x]])^2),{x,-50,0},PlotRange->All]

The grid and the box

We shall use a box [-L,L] and a grid composed by nP (internal) points. The unused points (Subscript[x, 0]) end (Subscript[x, nP+1]) coincide with -L and +L. The parameters nP and L will be passed to the functions below and can be changed. In a box of length Lbox the numbers of intervals is nP+1 and their width is Δx = Lbox/(nP+1).

In a finite box is better to smooth initial data in such a way that ψ vanishes exactly on the boundary. This can be done with multiplication by a function which is practically 1 almost everywhere and vanishes on the boundary as

smoothStep[x_,L_]:=-Tanh[x-L] Tanh[x+L];
Plot[smoothStep[x,10],{x,-10,10},PlotRange->{{-10,10},{0,1}}]

Time evolution: spatial grid

On a finite grid the evolution operator is approximated by

dU ≃ ((1 - i (1/2) H Δt)/(1 + i (1/2) H Δt));      group properties ->         ψ[n Δt] = (((dU))^n)ψ[0];

The potential on the grid is given by V[(Subscript[x, k])] the kinetic part of H (the laplacian operator) is approximated by finite difference approximation to second derivative, see notebook [NB-3.1] All this stuff is done by the function getdU already introduced in notebook [NB-3.4]

getdU[Δt_,V_,nP_,L_]:= Module[
{D2,indDiag,matV,h,x,Ham,U2,U1,U},
(* definition of the grid *)
h = N[2.L/(nP+1)]; 
x = - L + 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 *)
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}]

Let us choose our box, the number of points in the grid and the small time step Δt to be used to compute infinitesimal time evolution

L = 50; nP = 1000; Δt = 0.1;

The function getdU gives gridpoints, xpt,  and infinitesimal evolution operator, dU.

{xpt,dU} = getdU[Δt,V,nP,L];//Timing

Let us now produce an animation for the evolution of |ψ| . We choose numFrames as the total number of graphs to be produced and dtStepsPerFrame, as the number of infinitesimal Δt evolutions to be done between frames. These parameters depend on the problem at hand and the reader can change them.

numFrames = 75;
dtStepsPerFrame = 4;
mycolor = RGBColor[0.76817, 0.762127, 0.991638];

We define our initial wave function

ψGrid = ψ0[xpt]smoothStep[xpt,L];
maxψ = Max[Abs[ψGrid]];

and produce its time development. We store in matrix matPsi the value of ψ at each frame.

matPsi = Table[Do[ψGrid = dU.ψGrid;,{dtStepsPerFrame}];ψGrid,{j,1,numFrames}];//Timing

evolutionPlot = Table[ListPlot[Transpose[{xpt,Abs[matPsi[[j]]]^2}],Filling-> Axis,
Joined-> True,PlotRange-> {{-L,L},{0,1.1((maxψ)^2)}},Epilog-> {RGBColor[1,0,0],Rectangle[{aL,0},{aR,1}]}],
{j,1,numFrames}];

ListAnimate[evolutionPlot,AnimationRunning->False]

From the previous animation it is evident the reflection of the wave on the potential barrier. It is also evident a reflection on the (left) boundary, a spurious finite size effect which can be minimized allowing larger values of L. In this scale the transmitted wave is invisible, so let us select this part of ψ using a filter, i.e. a vector, gate, which put to 0 the wave function on the left, leaving only ψ after the barrier.

gate = N[UnitStep[xpt - aR]];
normaT =2Max[Abs[matPsi[[numFrames]]]*gate]^2;

The animation of the transmitted wave. The wave appears with some delay, due, roughly speaking, to the evolution from initial position x0P to the barrier. This graph is useful for energies lower than that used in this example, when transmitted wave can be almost invisible in previous animation.

ListAnimate[Table[Column[{Show[evolutionPlot[[j]],PlotRange-> {{-L,L},{0,normaT}}],j}],{j,1,numFrames}],AnimationRunning->False]

In this animation is even more evident the reflection on the right boundary.

Momentum analysis

In this section we show in detail how transmission coefficient works on a wave packet. Let us first of all choose a good transmitted wave, i.e. not effected by boundary reflection. A possible choice is the frame 62. This value obviously depends on the parameters used, in the general case the reader must look at the previous animations and make a suitable choice.

frameTransm =62;

ListPlot[Transpose[{xpt,((Abs[matPsi[[frameTransm]] gate])^2)}],Joined->True,PlotRange->{{0,L},{0,normaT}}]

trasmittedWave = matPsi[[frameTransm]]*gate;

The transmission coefficient is usually computed for plane waves and it depends on k, the wave number. For an incoming wave of the general form

(Subscript[ψ, in])[x] = (∫_(-∞)^∞) c[k]Exp[i k x] (dk/(2π))

Linearity of Schrödinger equation gives for the transmitted wave

(Subscript[ψ, Trasm])[x] = (∫_(-∞)^∞)(Subscript[a, T])[k] c[k]Exp[i k x] (dk/(2π))

The transmission coefficient,  T = | (Subscript[a, T])((|)^2) is directly related to the spectrum of the transmitted wave, i.e. the square of Fourier components

(Subscript[P, in])[k] = ((c[k])^2);    (Subscript[P, Trasm])[k] = |(Subscript[a, T])((|)^2)((c[k])^2)= T[k]* (Subscript[P, in])[k].

Fourier decomposition

There are several possibilities to extract the Fourier spectrum of a function. On a box of length Lbox with zero boundary conditions the permitted momenta are

(Subscript[k, j])= (π/Lbox)*j

One way to obtain the c[k]'s is to consider the Fourier series for a function extended in the interval of length 2Lbox in an odd way. In the doubled box the momenta generated by Fourier transform are multiple of 2π/(2 Lbox) = π/Lbox, the same as before.

The extension of ψ is easily constructed, is sufficient to add a zero value for the missing boundary point ( x=L, we remember that the grid contains only internal points) and join the original values of ψ in reverse order, changing its sign. The Mathematica function Fourier gives the Fourier coefficients. We normalize the spectrum to 1 so we do not care about an overall normalization factor.

cFFT = Fourier[Join[trasmittedWave,{0},-Reverse[trasmittedWave]]];

The allowed momenta are ( remember that Lbox = 2 L )

k = N[(π/(2L))Range[0,nP-1]];

We can write the spectrum, normalizing the result to its maximum for convenience.

spectrum=Take[((Abs[cFFT])^2),nP];
spectrum=(spectrum/(Max[spectrum]));
plotSpectrum=ListPlot[Transpose[{k,spectrum}],PlotRange->{{0,1.5` ϵ0},All},AxesLabel->{"k","P"}]

We can compare this result with theoretical expectations. In the text the transmission coefficient for a rectangular barrier of height (ϵ^2) /2 and width d has been computed

𝒯[k_,d_,ϵ_]:= Which[k≤ ϵ,((8 (k^2) ((k^2)-(ϵ^2)))/(8 (k^4)-8 (k^2) (ϵ^2)+(ϵ^4)-(ϵ^4) Cosh[2 d (Sqrt[-(k^2)+(ϵ^2)])])),
k>ϵ,((8 (k^2) ((k^2)-(ϵ^2)))/(8 (k^4)-8 (k^2) (ϵ^2)+(ϵ^4)-(ϵ^4) Cos[2 d (Sqrt[(k^2)-(ϵ^2)])]))]

The Fourier transform of the initial gaussian wave packet, is again a gaussian wave packet centered at k= k0

Integrate[Exp[I q0 x- α (((x-x0))^2)] Exp[-I q x],{x,-∞,∞},GenerateConditions-> False]

Then apart from an overall normalization factor we have

Pout = 𝒯[#1,aR-aL,ϵ0]*Exp[-(1/(2 (dk0^2)))(((#1 - k0))^2)]&/@ k;
Pout = Pout/Max[Pout];

The agreement with numerical result is very good. The red line indicates momentum k0, showing that transmitted wave is centered at higher momentum than incident packet.

Show[ListPlot[Transpose[{k,Pout}],Joined->True,PlotRange->{{0,1.5` ϵ0},All},AxesLabel->{"k","P"},Epilog->{Dashing[{0.02`,0.02`}],RGBColor[1,0,0],Line[{{k0,0},{k0,1}}]}],plotSpectrum]

Quantitatively we see that the error in the interesting region (where P ∼ 1) is of the order of  1%.

ListPlot[Transpose[{k,Pout-spectrum}],PlotRange->{{0,1.5` ϵ0},All}]

Time evolution: NDSolve

As shown in notebook [NB-3.4] time dependent  Schrödinger equation can be handled with Mathematica function NDSolve. Let us consider again zero boundary conditions in a box [-L,L] and define the form of the equation with L as parameter:

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

To compare with the previous computation we consider the same time interval:

tmax = Δt*numFrames*dtStepsPerFrame;

The solution of the equation and the production of an animation follows:

ψ=u/.First[NDSolve[equation[L],u,{x,-L,L},{t,0,tmax},PrecisionGoal-> 7]];//Timing

Animate[Plot[((Abs[ψ[x,t]])^2),{x,-L,L},Filling-> Axis,
PlotRange->{{-L,L},{0,1.1`}},Epilog->{RGBColor[1,0,0],Rectangle[{aL,0},{aR,1}]}],{t,0,tmax},AnimationRunning->False]

The qualitative behavior is as expected. In the following section we compare the two approach and discuss the precision of the results.

Comparison between NDSolve and finite grid

First of all let us compare the results for transmitted wave. We can construct the spectrum as before, taking ψ at the same t

ψNDS = ψ[xpt,frameTransm*dtStepsPerFrame*Δt]*gate;
cFFTNDS = Fourier[Join[ψNDS,{0},-Reverse[ψNDS]]];

spectrumNDS = Take[((Abs[cFFTNDS])^2),nP];
spectrumNDS = spectrumNDS/Max[spectrumNDS];

Show[ListPlot[Transpose[{k,Pout}],Joined->True,PlotRange->{{0,1.5` ϵ0},All},AxesLabel->{"k","P"},Epilog->{Dashing[{0.02`,0.02`}],RGBColor[1,0,0],Line[{{k0,0},{k0,1}}]}],ListPlot[Transpose[{k,spectrumNDS}],PlotRange->{{0,1.5` ϵ0},All},AxesLabel->{"k","P"}]]

ListPlot[Transpose[{k,Pout-spectrumNDS}],PlotRange->{{0,1.5` ϵ0},All}]

The precision is good as in the previous case.

Let us now consider the whole result.

ListPlot[Table[vec=ψ[xpt,j dtStepsPerFrame Δt];((Max[Abs[matPsi[[j]]-vec]])/(Max[Abs[vec]])),{j,1,numFrames}]]

We see that discrepancy grows with time and becomes bigger when the wave hits the barrier. Due the coincidence of the results for transmitted wave the discrepancy is due to the complicated reflection pattern. This can be shown for example plotting the difference in the real part of ψ.

Show[GraphicsGrid[Partition[Table[ListPlot[Transpose[{xpt,Re[matPsi[[j]]]-Re[ψ[xpt,j dtStepsPerFrame Δt]]}],PlotRange->All,PlotLabel->"j="<>ToString[j]],{j,Floor[(numFrames/4)],numFrames,Floor[(numFrames/4)]}],2]],ImageSize -> 72 8]

The question arises: which method is more precise? As expected Mathematica internal function NDSolve is better: its internal algorithm is based on a method similar to the one used on the grid but allows a change in running time for time steps Δt, in this way the solution can follows the rapid variation at the reflection on the barrier. To prove this point let us redo the finite grid computation with a smaller Δt and a greater dtStepsPerFrame:

scale = 2;
{xpt,dU} = getdU[Δt/scale,V,nP,L];
dtStepsPerFrame=2*scale;

ψGrid = ψ0[xpt]smoothStep[xpt,L];
matPsiNew ={};
(Table[Do[ψGrid = dU.ψGrid;,{dtStepsPerFrame}];
matPsiNew = Append[matPsiNew,ψGrid];
,
{j,1,numFrames}];)

ListPlot[Table[vec=ψ[xpt,((j dtStepsPerFrame Δt)/scale)];((Max[Abs[matPsiNew[[j]]-vec]])/(Max[Abs[vec]])),{j,1,numFrames}]]

The discrepancy has been indeed lowered. Even better can be done with more grids points ( nP ). The check is left to reader

In conclusion, as expected, Mathematica procedure NDSolve is to be preferred.

Exercises

Use the function NDSolve with periodic boundary conditions to solve the problem exposed in this notebook.

Compute < x > and < p > for transmitted and reflected waves.