Quantum Dynamic time evolution using the Split-Operator Fourier Transform (SOFT) algorithm

Some ideas what you can do with the code in this notebook:

- Compare the FFT implemented by mathematica and self-written summation version (which is better?)
-Compare different split-ups (Theory tells us, that the symmetric Trotter splitup presented in the video should perform better than other splitups). I started implementing another, presumably more “stupid” splitup. You can use this to compare different splitups with the symmetric one
-Change the grid-spacing, change the time steps, what changes and why?
-Change the grid for p. Does the code still work? Why? Why not? 
-Change the potential (already defined are a harmonic and a wall potential - change them or try different ones)
-Change the form of the wavepackets or use two of them (propagate the sum of two such Gaussians) what happens? 
....

The Problem

The approach - setup

First, we have to define our grids. We define lower and upper limits and resolutions for the grid in the coordinate and the grid in the momentum space.

xmin = -5;
xmax = 5;
steps =  2000;
m = 0.5;

gridx = Array[#&,steps,{xmin,xmax}];
gridp = Array[#&,steps,{xmin,xmax}];

For the Fourier transform steps, we need to know the spacing of two neighboring elements in the list:

spacingx =Abs@(gridx[[1]]-gridx[[2]]) ;(*dx*)

spacingp =Abs@(gridp[[1]]-gridp[[2]]) ; (*dp*)

For the time evolution, we have to define the width of our time steps as well as the runtime.

tmax = 100;
dt = 0.001;

Fourier transforms

In any case - i.e. independent of the splitting approach - we have to carry out two kinds of Fourier transformations, namely from the x to the p space and of course the other way round.  To keep it as close as possible to the formula, we use a Table[] (although mathematica has a function for DFT (tutorial/FourierTransforms).

fourierxp[phix_,dx_] := Table[
Total[phix*Exp[-I*p*gridx]]*dx,{p,gridp}
]

fourierpx[phip_,dp_] := Table[
Total[phip*Exp[I*x*gridp]]*dp,{x,gridx}
]/(2*Pi)

Symmetric trotter expansion

symmTrotEv[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := 
Module[
				{psi = WF,
				gx = gridx,
				gp = gridp,
				dx = spacingx,
				dp = spacingp,
				dt = spacingt,
				pot = potential
				},
				psi = Exp[-I*pot*dt/2] * psi; (*in x space*) 
				psi = fourierxp[psi,dx]; 
				psi = Exp[-I*gp^2 *dt/(2*m)]*psi; (*in p space*)
				psi = fourierpx[psi,dp];
				psi = Exp[-I*pot*dt/2]*psi (*in x space and output *)
			];

symmTrotEv2[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := 
Module[
				{psi = WF,
				gx = gridx,
				gp = gridp,
				dx = spacingx,
				dp = spacingp,
				dt = spacingt,
				pot = potential
				},
				psi = Exp[-I*pot*dt/2] * psi; (*in x space*) 
				psi = Fourier[psi]; 
				psi = Exp[-I*gp^2 *dt/(2*m)]*psi; (*in p space*)
				psi = InverseFourier[psi];
				psi = Exp[-I*pot*dt/2]*psi (*in x space and output *)
			];

symmTrotEv2[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := 
Module[
				{psi = WF,
				gx = gridx,
				gp = gridp,
				dx = spacingx,
				dp = spacingp,
				dt = spacingt,
				pot = potential
				},
				psi = Exp[-I*pot*dt/2] * psi; (*in x space*) 
				psi = Fourier[psi]; 
				psi = Exp[-I*gp^2 *dt/(2*m)]*psi; (*in p space*)
				psi = InverseFourier[psi,dp];
				psi = Exp[-I*pot*dt/2]*psi (*in x space and output *)
			];

“Stupid” trotter expansion

stTrotEv[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := 
Module[
				{psi = WF,
				gx = gridx,
				gp = gridp,
				dx = spacingx,
				dp = spacingp,
				dt = spacingt,
				pot = potential
				}, 
				psi = Exp[-I*pot*dt]*psi (*in x space and output *);
				 psi = fourierxp[psi,dx]; 
				psi = Exp[-I*gp^2 *dt/(2*m)]*psi; (*in p space*)
				psi = fourierpx[psi,dp]	
			];

stTrotEv2[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := 
Module[
				{psi = WF,
				gx = gridx,
				gp = gridp,
				dx = spacingx,
				dp = spacingp,
				dt = spacingt,
				pot = potential
				}, 
				psi = Exp[-I*pot*dt]*psi (*in x space and output *);
				 psi = Fourier[psi]; 
				psi = Exp[-I*gp^2 *dt/(2*m)]*psi; (*in p space*)
				psi = InverseFourier[psi]	
			];

Potential

harmonicpotential[x_] = 0.05*x^2;
harmonicpotgrid = harmonicpotential[#]&/@gridx;

wallpoential[x_] :=Piecewise[ {{1000000,-3.2≤ x≤ -2.8},{0,x<-3.2},{0,x>-2.8}}]

wallpoentialgrid = wallpoential[#]&/@gridx;

Plot[wallpoential[x] ,{x,-5,10}]

Gaussian wave packet

To define a wave packet, we have to specify where it is centered and how steep it is.

xcenter = 0; (*must be larger than xmin!*)
gaussparam = 0.5;

wp = Exp[-(gridx-xcenter)^2/(2*gaussparam^2)]*(1/Sqrt[gaussparam*Sqrt[Pi]]);

norm = Sqrt[Total[Conjugate[wp]*wp]]

wp= wp/norm;

ListPlot[wp,PlotRange->Full]

Total[Conjugate[wp]*wp]

Convergence

Two often used convergence criteria are normalization and energy, if the time evolution is “correct” both should remain conserved.

normalization[f_,t_]  := Module[{wpn=f,time=t,norm},
norm= Total[Conjugate[wpn]*wpn];
{time,norm}
]

energy[f_,t_]  := Module[{wpe=f,time=t,ekin, epot,etot,k},
epot= Total[Conjugate[wpe]*wallpoentialgrid*wpe];
k= 2*Pi*spacingp;
ekin = Total[Conjugate[wpe]*0.5/m*k^2*wpe];
etot =Re@(ekin+epot);
{time,etot}
]

normalization[wp2,0]

Time evolution

Now, we run all the functions which we have written to study the time evolution.

t= 0; (*Initialization*)

(*symmTrotEv[WF_,gridx_,gridp_,spacingx_,spacingp_,spacingt_,potential_] := *)

convergence = {};

energyconv= {};

AbsoluteTiming@While[t <tmax,
wpint = symmTrotEv2[wp,gridx,gridp,spacingx,spacingp,dt,wallpoentialgrid];
wp = wpint;
AppendTo[convergence,normalization[wp,t]];
AppendTo[energyconv,energy[wp,t]];
t=t+dt;
]

convergencetinystep =energyconv;

convergencemidstep =energyconv

movie = Dynamic[
GraphicsRow[{
Show[
ListPlot[Partition[Riffle[gridx,Re@wp],2],Joined->True,PlotRange->{-0.10,0.15},AxesLabel-> {x,"ψ / |ψ((|)^2) / V"}],
ListPlot[Partition[Riffle[gridx,norm*Conjugate[wp]*wp],2],Joined->True,PlotRange->{-0.05,0.15},PlotStyle->Red],
ListPlot[Partition[Riffle[gridx,wallpoentialgrid],2],Joined->True,PlotRange->{-1,20}] ],
ListPlot[convergence,AxesLabel->{"time","norm"}],ListPlot[energyconv,PlotRange->All,AxesLabel->{"time","energy"}]},ImageSize->Full
]]

Stupid time develop

wp2= Exp[-(gridx-xcenter)^2/(2*gaussparam^2)]*(1/Sqrt[gaussparam*Sqrt[Pi]]);

norm2= Sqrt[Total[Conjugate[wp2]*wp2]]

wp2 =  wp2/norm2;

t2= 0;

convergence2 = {};

energyconv2= {};

AbsoluteTiming[While[t2 <tmax,
wpint2 = stTrotEv2[wp2,gridx,gridp,spacingx,spacingp,dt,wallpoentialgrid];
wp2 = wpint2;
AppendTo[convergence2,normalization[wp2,t2]];
AppendTo[energyconv2,energy[wp2,t2]];
t2=t2+dt;
]]

Dynamic[t2]

Dynamic[
GraphicsRow[{
Show[
ListPlot[Partition[Riffle[gridx,Re@wp2],2],Joined->True,PlotRange->{-0.05,0.15}],
ListPlot[Partition[Riffle[gridx,norm*Conjugate[wp2]*wp2],2],Joined->True,PlotRange->{-0.05,0.15},PlotStyle->Red],ListPlot[Partition[Riffle[gridx,wallpoentialgrid],2],Joined->True,PlotRange->{-1,10}] ]
,ListPlot[convergence2],ListPlot[energyconv2,PlotRange->Full]},ImageSize->Full
]

]

Show[ListPlot[convergence2,PlotStyle->Red],ListPlot[convergence],ImageSize->Large]

Show[ListPlot[convergencemidstep,AxesLabel->{"time","energy"}],ListPlot[energyconv,PlotStyle->Red,PlotRange->Full],ImageSize->Large]