Time evolution

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the time evolution of a wave packet in an external potential by three methods

In particular study the tunneling  process for a double well potential.

Utilities

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

We shall work on a finite space domain, [-L,L]. It is often useful to smoothly join to zero a function at the boundaries, this can be done with the following function:

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

a second possibility is to use Mathematica function UnitStep: UnitStep[L-x,x+L];

General remarks

In quantum mechanics time evolution is described by the time-dependent Schrödinger equation:

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

Units

It is convenient and important to choose dimensionless units in problem solving. This allow to understand which parameters are relevant. A system of units must fix 3 scales (as in CGS system). In our problems two of them are natural, we put m= 1, and ℏ=1 for units of mass and action. The third one typically is chosen as a characteristic length (or time) of the problem. We will mainly study evolution of gaussian packets

ψ(x) = (1/((π^(1/4))(Sqrt[δx])))Exp[-(( (x^2))/(2 ((δx)^2)))];   width:  ((Δx)^2) = ((δx)^2)/2;

It is natural to put as unit of length δx. With such a choice for example we have for derived physical quantities:

energy: ((ℏ^2)/(m ((δx)^2)));      time:(( m ((δx)^2))/ℏ);  etc.

Analytical results for the free particle

In simple cases ( quadratic Hamiltonians ) it is possible to write an analytical expression for the fundamental solution (Green's function) of the Schrödinger equation and reduce to a quadrature the evolution of a given state. This procedure will be explored later in this course. For the elementary methods exposed in this notebook it is sufficient to remember some trivial results for the free particle evolution. If an initial state is given as superposition of plane waves:

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

Its time evolution is given by

ψ[x,t] = (∫_(-∞)^∞)f[k]Exp[i k x]Exp[-i ((ℏ (k^2))/(2m))t]((dk)/(2π))

In particular a normalized gaussian packet moving with mean momentum q is described by

ψ[x,0]=(1/((π^(1/4))(Sqrt[δx])))Exp[-(( (x^2))/(2 ((δx)^2)))]Exp[i q x];      f[k]=  (π^(1/4))(Sqrt[2 δx]) (e^(-(1/2) (((k-q))^2) ((δx)^2))) .

Performing the k integration in () one gets

ψ[x,t]=((e^((-(q^2) t+2 q x+i (x^2)-2 (q+i x) (Subscript[x, 0])+i (Subscript[x, 0]^2))/(2 (-i+t))))/((π^(1/4)) (Sqrt[1+i t])))

The reader is strongly encouraged to check above formulas, maybe using Mathematica .

The probability distribution, which will be used as a check in the following, is given by (in natural units):

densityFreeParticle[x_,t_,q_,x0_]:=((E^(-((((q t-x+x0))^2)/(1.+(t^2)))))/((Sqrt[π]) (Sqrt[1.+(t^2)])))

Numerical procedures

We give 3 different albeit related methods for studying time evolution:

Details on each method will be given below. The first two methods will be checked for the free particle case. The methods 1 and 3 will be applied to a double well potential. All the procedures will work for a large class of potentials, which method is better depends on the problem.

Discretization method

In method 1) we consider a time discretization in small intervals Δt. The multiplication property of the exponential (i.e. group properties of time translation) implies

t= N Δt  -> U = (((dU))^N);     dU = Exp[- i H Δt];

For small intervals Δt we can approximate the exponential by:

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

In numerical analysis the approximation () is known as the Cranck-Nicholson scheme. From our point of view the main point is that approximation () is explicitly unitary,i.e.  (dU^(+ ))= (dU^(-1)) as the reader can easily check.

At fixed time we discretize space with a grid of nP points, as in notebook [1]. In this approximation H becomes a matrix and eq.() provides the solution via matrix multiplication. See notebook [1] and below for details.

Space grid

In a numerical method it is necessary to work with finite intervals:  we shall take  x ∈ [-L,L]. This interval will be divided in nP+1 intervals, with nP internal points. The boundary condition ψ[±L] = 0 is imposed and only ψ values at internal points are to be found. If (Subscript[x, k])= k h, k=1... nP denote the points in the grid, the potential is discretized by a diagonal matrix constructed withV[[k]] = V[ (Subscript[x, k])]. Kinetic term comes from a discrete version of second derivative:

f''[x[j]] ≃ ( f[j+1] - 2 f[j] + f[j-1])/(h^2);

An example for nP = 5, taking account the boundaries conditions, is:

ψ''≃(1/(h^2)) D2.ψ ;                         D2 ψ= (({{-2, 1, 0, 0, 0}, {1, -2, 1, 0, 0}, {0, 1, -2, 1, 0}, {0, 0, 1, -2, 1}, {0, 0, 0, 1, -2}})) (({{ψ[[1]]}, {ψ[[2]]}, {ψ[[3]]}, {ψ[[4]]}, {ψ[[5]]}}));  (nP = 5);

Time grid

We will study a finite time  interval [0,tmax]. This interval will be divided in a number of numFrames small intervals, at the end of each interval ψ will be computed and eventually plotted. Each small interval will be taken of the size δt = Δt * dtStepsPerFrame, i.e. a number of dtStepsPerFrame of relation () are performed, each with a time-step Δt, for each small time - interval. This separation of Δt and dtStepsPerFrame is done to keep separated the approximation in ()  from the number of actual frames to be taken (for example to do a plot).

The evolution function dU

All the previous work is done by the following function dU. The function dU[Δt,V,nP,L] needs as input the time step, the potential function, the number of grid points ( internal to the interval) and the size of space-grid. Here we put the system in a box [-L, L], the grid contains nP+1 intervals and h = 2L/(nP+1).  The output gives the list of grid points ( x ) and the one step time evolution operator, dU. The If instruction handle the case of free particles, for which V[x] = 0. For the function SparseArray the reader is referred to Mathematica Help.

getdU[Δt_,V_,nP_,L_]:= Module[
{D2,indDiag,matV,h,x,Ham,dU,U2,U1},
(* 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 *)
If[NumberQ[V],matV = 0,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}]

Time evolution

We consider the case of a free particle, with the initial wave function (in natural units):

ψ0[x_]:= (1/(π^(1/4)))Exp[-(( (x^2))/(2 ))]Exp[I  x];

Let us define the space and time grids

L = 50; nP = 1000; 
Δt = 0.1; dtStepsPerFrame = 2; tmax=15;
numFrames = Floor[tmax/(dtStepsPerFrame*Δt)];

The following steps will be repeatedly used below and in others notebooks.

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

The first instruction define the grid and the infinitesimal time evolution operator dU.

The second smoothens the initial function in our interval (to get exactly zero boundary conditions).  maxψ is computed for use in a plot and finally the relation ψ[t+δt] = U[δt] ψ[t] is used recursively.

The result is a matrix, matPsi, In each row of the matrix, matPsi[[j]], the values of ψ at the points of the grid are stored.

Let us make an animation for the evolution of particle density of the packet. For clarity the plot of initial distribution is done and Prepended to the list of evolutionPlot.

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

We can check our result against the analytical function densityFreeParticle.

ListPlot[Table[vec=densityFreeParticle[xpt,j dtStepsPerFrame Δt,1,0];((Max[((Abs[matPsi[[j]]])^2)-vec])/(Max[Abs[vec]])),{j,1,numFrames}],PlotRange->All,AxesLabel->{"frames","relative error"}]

As we see there is an error of about 10 per cent at long times. The results can be improved by taking a smaller Δt and a greater number of points, with a cost in memory and time. As an example let us take nP = 2000.

Some variation of the parameters

The computations of this subsection take some more time and can be skipped in a first reading.

Let us store our old result for later use

matPsiOld = matPsi;

nP = 2000;

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

ListPlot[Table[vec=densityFreeParticle[xpt,j dtStepsPerFrame Δt,1,0];((Max[((Abs[matPsi[[j]]])^2)-vec])/(Max[Abs[vec]])),{j,1,numFrames}],PlotRange->All,AxesLabel->{"frames","relative error"}]

We can alternatively take half time intervals, doubling dtStepsPerFrame to take always the same number of output values for ψ.

nP = 1000; 
Δt = Δt/2; dtStepsPerFrame = 2*dtStepsPerFrame;

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

ListPlot[Table[vec=densityFreeParticle[xpt,j dtStepsPerFrame Δt,1,0];((Max[((Abs[matPsi[[j]]])^2)-vec])/(Max[Abs[vec]])),{j,1,numFrames}],PlotRange->All,AxesLabel->{"frames","relative error"}]

Let us restore our old values.

matPsi = matPsiOld;
nP=1000;Δt = 0.1; dtStepsPerFrame = 2;
Clear[matPsiOld];

NDSolve

The use of Mathematica function NDSolve is the most direct approach to the study of time evolution. Internally Mathematica uses (in most cases)  a method similar to a space - time grid, except that more sophisticated methods are employed to solve time dependent evolution at each space grid point. To use the full power of NDSolve is not so simple for a non expert user, and at least in this notebook we try only the simplest possibilities. Obviously the reader is strongly invited to consult the ample documentation in Mathematica Help, especially the internal link Advanced Documentation.

The first observation is on the smoothing procedure. Mathematica tries to do calculations as best as possible: smoothing a function as previously done in regions where its value is very small, typically (10^(-40)) is  time consuming and meaningless from the side of the precision of the methods used. It is better to just cut the function with smoothStep:

ψ0ND[x_]:= (1/(π^(1/4)))Exp[-(( (x^2))/(2. ))]Exp[I  x]smoothStep[x,L]

We have not taken this attitude in the previous method because for algebraic calculations previous smoothing is more stable.

With zero boundary conditions the solution of our problem is just two lines code. First define an equation (we leave L as parameter if the reader wants to change L (in this case can be necessary to change the smoothing).

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

The solution is given by:

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

The reader has certainly noticed that the procedure is very fast with respect to the previous code. Let us analyze the error at the same grid points

error[ψ_]:=ListPlot[Table[vec=densityFreeParticle[xpt,j dtStepsPerFrame Δt,1.`,0.`];vec2=((Abs[ψ[xpt,j dtStepsPerFrame Δt]])^2);((Max[vec2-vec])/(Max[Abs[vec]])),{j,1,numFrames}],PlotRange->All]

error[ψ]

We see that the error is smaller than the previous case, but grows around 60 frames. This means that for longer times this simple method must be adapted, working with the options of NDSolve.

An even faster method is to change the boundary conditions, putting periodic boundary conditions in NDSolve.

equationPeriodic[L_]:= {I(∂_t)u[x,t]==-1/2(∂_(x,x))u[x,t] ,u[x,0]==ψ0ND[x],u[L,t]==  u[-L,t]}

ψ=u/.First[NDSolve[equationPeriodic[L],u,{x,-L,L},{t,0.,tmax},
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", DifferenceOrder->"Pseudospectral" }}]];//Timing

The situation is better:

error[ψ]

This method is quite stable in the present case and a greater precision can be required:

ψ=u/.First[NDSolve[equationPeriodic[L],u,{x,-L,L},{t,0.,tmax},PrecisionGoal-> 10,AccuracyGoal-> 12,
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", DifferenceOrder->"Pseudospectral" }}]];//Timing

error[ψ]

The option  DifferenceOrder->"Pseudospectral" refers to a method analogous to Fourier transform, the reader is referred to Mathematica help for more details.

Harmonic oscillator

Let us now repeat the procedure for the harmonic oscillator potential.

L=10;
tmax = 2.*2π;

V[x_]:= (1/2)(x^2);

We define the initial state as a gaussian packet, with mean momentum 1 (in natural units) and dispersion Δx2.

Δx2=0.4;
ψ0ND[x_]:= (1/(π^(1/4)))Exp[-(( (x^2))/(4.Δx2 ))]Exp[I  x]smoothStep[x,L]

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

ψ=u/.First[NDSolve[equationPeriodic[L],u,{x,-L,L},{t,0.,tmax},
PrecisionGoal-> 10,
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", DifferenceOrder->"Pseudospectral" }}]];

Δt=0.1;
times = Range[0,tmax,Δt];

Here is an animation for the evolution of the state. We can appreciate the oscillation in the width of the state.

xgridPlot = Range[-4,4,0.1];
ListAnimate[ Table[
ListPlot[Transpose[{xgridPlot,Abs[ψ[xgridPlot,t1]]^2}],Filling-> Axis,
Joined-> True,PlotRange-> {{-4,4},{0,0.6}}],{t1,times}],AnimationRunning->False]

Coherent state

We know (see text) that in natural units a dispersion Δx2 = 1/2 would correspond to a coherent state. Let us verify the difference in the evolution

Δx2=1./2;

ψ0ND[x_]:= (1/(π^(1/4)))Exp[-(( (x^2))/(4.Δx2 ))]Exp[I  x]smoothStep[x,L]

ψ=u/.First[NDSolve[equationPeriodic[L],u,{x,-L,L},{t,0.,tmax},
PrecisionGoal-> 10,
Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", DifferenceOrder->"Pseudospectral" }}]];

xgridPlot = Range[-4,4,0.1];
ListAnimate[ Table[
ListPlot[Transpose[{xgridPlot,Abs[ψ[xgridPlot,t1]]^2}],Filling-> Axis,
Joined-> True,PlotRange-> {{-4,4},{0,0.6}}],{t1,times}],AnimationRunning->False]

Eigenfunction expansion

From quantum mechanical point of view the problem of time evolution is reduced to the diagonalization of the Hamiltonian. If  ψB[k,x] are the eigenfunctions, E[k] the energy eigenvalues and ψ0[x] the initial state we have directly (we put as usual ℏ = 1):

ψ0[x] = ∑ c[j] ψB[j,x];       ψ[x,t] = ∑ c[j] Exp[-i E[j] t ]ψB[j,x]

We already know different methods to obtain approximate eigenstates, so we can find an approximate solution by a truncation of previous equation. We shall study quite sophisticated methods to obtain E[j] and ψB[j,x], in this notebook we use the very simple procedure developed in notebook [3]. Let us recall the main points

ψEigenstate[[k,x]]= (∑_(j=1)^dim) eigV[[k,j]]ψH[j,x]

c[[k]] = ∫ ψ0[x]*ψEigenstate[[k,x]] =  (∑_(j=1)^dim) eigV[[k,j]]∫ ψ0[x]*ψH[j,x]

A model study

To take a nontrivial case let us consider a double well potential

V[x_]:= (1/2)g ((((x^2)-(μ^2)))^2);

The potential has two local minima. Quantum mechanically we expect that a particle confined in one well for t=0 will tunnel and start to oscillate between the two miima. Expanding around x=μ, for x = μ+ξ the potential has a quadratic part 1/2 4g (μ^2)(ξ^2) i.e. e frequency ω = 2μ(Sqrt[g]). An approximate ground state is given by the following function (with c=1). (A free parameter c is given to sharpen the initial state if necessary).

ψ0[x_]:= (((c g))^(1/8))(((2 μ/π))^(1/4))Exp[-(Sqrt[c g])μ (((x-μ))^2)];

Let us choose the parameters g and μ

μ=1.3; g=1.; c= 1;

The state is normalized to unity:

NIntegrate[((ψ0[x])^2),{x,-∞, ∞}]

We note that for these parameters the mean energy is lower than local maximum of the potential g (μ^4)/2 , i.e. transition to the left hand side of the potential is classically forbidden.

{meanEnergy=NIntegrate[(1/2)((((∂_x)ψ0[x]))^2)+ ((ψ0[x])^2)(1/2)g ((((x^2)-(μ^2)))^2),{x,-∞, ∞}],g (μ^4)/2}

Numerical time evolution

Hamiltonian and 0-order eigenstates

The function XPmatrices has been defined in notebook [3]

ℏ =1; m=1; ω = 1;

XPmatrices[dim_,parity_Integer:-1]:=
Module[{x2array,p2array,x4array,x2,p2,x4,x1,x1array,p1,p1array}, 
Clear[i,n];
singlepower= False;
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;singlepower = True];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[singlepower == True,
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
p1= SparseArray[{{i_,k_}/;k-i==1->-I(Sqrt[(n+1)/2])},{dim,dim}];
p1 = (p1-Transpose[p1]);
aOp = SparseArray[{{i_,k_}/;k-i==1->(Sqrt[n+1])},{dim,dim}];
adagOp = Transpose[aOp];];
If[parity == 0 || parity == 1,
Normal[#1]&/@{m ω ℏ p2,(ℏ/(m ω))x2,((((ℏ/(m ω))))^2)x4},
Normal[#1]&/@{(Sqrt[m ω ℏ])p1,m ω ℏ p2,(Sqrt[(ℏ/(m ω))])x1,(ℏ/(m ω))x2,((((ℏ/(m ω))))^2)x4,aOp,adagOp}]];

We choose a dimension (the number of independent bases functions)

dim = 50;

and compute the necessary operators:

{P,P2,Q,Q2,Q4,aOp,adagOp} = N[XPmatrices[dim]];

The Hamiltonian of the system is easily defined

H = P2/2  + g/2 (Q4 - 2 (μ^2)Q2 + (μ^4)IdentityMatrix[dim]);

while the harmonic oscillator basis on which H has been constructed are given by

ψH[n_,x_]:= (1/(π^(1/4)))(1./(Sqrt[(2^n)(n)!]))HermiteH[n,x]Exp[-x^2/2];

Diagonalization

{ene,eigV} = Eigensystem[H];
ord = Ordering[ene];
{ene,eigV} = {ene[[ord]],eigV[[ord]]};
eigV = eigV/(Norm/@eigV);

Off[NIntegrate::"ncvb"]

coeffHermite=Table[NIntegrate[ψ0[x]*ψH[n,x],{x,-∞,+∞}],{n,0,dim-1}];

On[NIntegrate::"ncvb"];

We write again the grid parameters in case the reader has made some changes in previous parts.

L = 5; nP = 1000; 
Δt = 0.1; dtStepsPerFrame = 2; tmax=40;
numFrames = Floor[tmax/(dtStepsPerFrame*Δt)];
h = N[2.L/(nP+1)]; 
xpt = - L + h Range[nP];

We select a subgrid for comparison with previous results and for plotting.

indexGrid=Range[1,Length[xpt],10];
xGrid = xpt[[indexGrid]];

We now compute harmonic oscillator functions on the grid and collect the values in matrix matPsiHermite.

matPsiHermite = Table[ψH[n,xGrid],{n,0,dim-1}];

A simple matrix multiplication gives eigenstates, see eq.(4.)

matPsiBasis = eigV.matPsiHermite;

and coefficient in this new base:

coeffBasis = eigV.coeffHermite;

As a check let use compare the values of original state, ψ0, with the one reconstructed from matPsiBasis

Max[Abs[coeffBasis.matPsiBasis - ψ0[xGrid]]]

Evolution

maxψ = Max[Abs[ψ0[xGrid]smoothStep[xGrid,L]]];
ListAnimate[ Table[
phasesEnergy = Exp[-I  ene*j*Δt*dtStepsPerFrame];
ψ = (phasesEnergy*coeffBasis).matPsiBasis;
ListPlot[Transpose[{xGrid,Abs[ψ]^2}],Filling-> Axis,Joined-> True,PlotRange-> {{-L,L},{0,1.3((maxψ)^2)}}],
{j,1,numFrames}],AnimationRunning->False]

The method works quite well. From the animation it is clear that two kind of oscillations are present: a short period one corresponding to the oscillation of the bump around the local minimum, a long period one connecting the two minima.

We leave as an exercise to the reader to make a comparison with the discretization technique.

Exercises

Consider a potential well V[x] = - (Subscript[U, 0]) for |x| < a and V[x] = 0 for |x| ≥ a. Take as the initial state a wave function confined in [-a,a] and study its evolution. What happens if <H> ≃ 0 and negative? Discuss the results

Apply the grid method to the double well potential and compare the results with eigenfunctions expansion.