time dependent Schröidinger Equation

Welcome back!

I prepare here an introduction to solve time-dependent Schröidinger equation and added my code at the end. If you are interested more professional resources following link may help
Wavepacket Dynamics
But less likely you may learn about the details of their code. So here is your chance to start playing one of the most exciting and instructive problem in Quantum Physics. Below you see a wavepacket scattering from a potential at the center of the graph.

Cheers!

During my research, I have to solve 2D Schroedinger equation both time dependent and independent for specific potentials. In following, I explain the idea behind how to solve time-dependent Schroidinger equation. It is a useful research tool and it is also instructive.
One of my research area in general can be categorized as transport in nanostructures. For that I want to solve time- dependent Schroidinger equation.

$i\mathrm{\hslash \; <!--\; planck\; constant\; over\; two\; pi\; -->}\frac{\partial \; <!--\; partial\; differential\; -->\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}{\partial \; <!--\; partial\; differential\; -->t}=-\frac{{\mathrm{\hslash \; <!--\; planck\; constant\; over\; two\; pi\; -->}}^2}{2m}\times \; <!--\; multiplication\; sign\; -->\left(\frac{{\partial \; <!--\; partial\; differential\; -->}^2\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}{\partial \; <!--\; partial\; differential\; -->x^2}+\frac{{\partial \; <!--\; partial\; differential\; -->}^2\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}{\partial \; <!--\; partial\; differential\; -->y^2}\right)+V\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}$

You can separate time dependence (separation of variable)

$\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}\left(x,y,t\right)=e^{-iHt/\mathrm{\hslash \; <!--\; planck\; constant\; over\; two\; pi\; -->}}\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}\left(x,y\right)H=-\frac{{\mathrm{\hslash \; <!--\; planck\; constant\; over\; two\; pi\; -->}}^2}{2m}\times \; <!--\; multiplication\; sign\; -->\left(\frac{\partial \; <!--\; partial\; differential\; -->}{\partial \; <!--\; partial\; differential\; -->x^2}+\frac{\partial \; <!--\; partial\; differential\; -->}{\partial \; <!--\; partial\; differential\; -->y^2}\right)$

For a computer simulation first I choose a unit system where energy is measured in units of something so that all constants are 1in Schroidinger eqn. This is usually the first step doing simulations and known as atomic units. In this unit wave function can be propagated simply

${\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}^{t+\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}=e^{-iH\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}^t$

For this problem it turns out to be the form that we write above is not stable numerically. we will use following form,

$e^{iHt/2}\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}=e^{-iHt/2}\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}\Rightarrow \; <!--\; rightwards\; double\; arrow\; -->\left(1+\frac{iHt}{2}\right)\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}\left(n,m,t\right)=\left(1-\frac{iHt}{2}\right)\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}\left(n,m,0\right)$

So we first multiply with half of the operator than we use ex pansion of the exponential and take the first term.
This is much stable numerically. The spatial part is also discritized as usual,

$\frac{{\partial \; <!--\; partial\; differential\; -->}^2\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}{\partial \; <!--\; partial\; differential\; -->x^2}=\frac{{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n+1}-2{\mathrm{\Psi \; <!--\; greek\; capital\; letter\; psi\; -->}}_n+{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n-1}}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}{x}^2},\frac{{\partial \; <!--\; partial\; differential\; -->}^2\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}{\partial \; <!--\; partial\; differential\; -->y^2}=\frac{{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{m+1}-2{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_m+{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{m-1}}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}{y}^2}$

Finally we get following eqn

$\begin{array}{}\left[\left(1-i{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_x\frac{\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}{2{\mathrm{\epsilon \; <!--\; greek\; small\; letter\; epsilon\; -->}}^2}\right)\times \; <!--\; multiplication\; sign\; -->\left(1-i{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_y\frac{\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}{2{\mathrm{\epsilon \; <!--\; greek\; small\; letter\; epsilon\; -->}}^2}\right)\times \; <!--\; multiplication\; sign\; -->e^{iVt/2}{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n,m}^{j+1}=\left(1+i{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_x\frac{\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}{2{\mathrm{\epsilon \; <!--\; greek\; small\; letter\; epsilon\; -->}}^2}\right)\times \; <!--\; multiplication\; sign\; -->\left(1+i\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}y\frac{\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}}{2{\mathrm{\epsilon \; <!--\; greek\; small\; letter\; epsilon\; -->}}^2}\right)\right]\times \; <!--\; multiplication\; sign\; -->e^{-iVt/2}{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n,m}^j\\ {\mathrm{△\; <!--\; latex:\; \backslash triangle\; -->}}_x{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n,m}^j={\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n+1,m}^j-2{\mathrm{\Psi \; <!--\; greek\; capital\; letter\; psi\; -->}}_{n,m}^j+{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n-1,m}^j\end{array}$

where we defined central derivative operator for x and y. Here ✓ is the static potential as above. delta and epsilon are time and space increments. This eqn has 3 indicies, it is useful to define a new variable and solve it in 2 parts.

$\left(1-i\mathrm{\lambda \; <!--\; greek\; small\; letter\; lambda\; -->}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_x\right){\mathrm{\varphi \; <!--\; greek\; phi\; symbol\; -->}}_n^j,m=\left(1+i\mathrm{\lambda \; <!--\; greek\; small\; letter\; lambda\; -->}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_y\right){\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n,m}^j$

Now,

$\left(1-i\mathrm{\lambda \; <!--\; greek\; small\; letter\; lambda\; -->}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_y\right)\times \; <!--\; multiplication\; sign\; -->e^{iV\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}/2}{\mathrm{\psi \; <!--\; greek\; small\; letter\; psi\; -->}}_{n,m}^{j+1}=\left(1+i\mathrm{\lambda \; <!--\; greek\; small\; letter\; lambda\; -->}{\mathrm{\Delta \; <!--\; greek\; capital\; letter\; delta\; -->}}_x\right)\times \; <!--\; multiplication\; sign\; -->e^{-iV\mathrm{\delta \; <!--\; greek\; small\; letter\; delta\; -->}/2}{\mathrm{\varphi \; <!--\; greek\; phi\; symbol\; -->}}_{n,n}^j$

first we use an initial condition and assume that at time step j we know the solution so phi can be calculated. So left hand side of above eqn is a known vector. Therefore rig a t hand side also calculated using 3 diagonal matrix form.

function M=tprop(J)
sigma=0.08;
a=0.08;
xv=0.4;
yv=0.4;
eps=0.01;
eta=0.5*eps^2;
k0=90;
kx=64;
ky=64;
Vmax=2*k0^2;
x=linspace(-0.5,1.5,J);
y=linspace(-0.5,1.5,J);
fi0=zeros(J,J);
V=zeros(J,J);
x0=0.25;
y0=0.25;
for ii=1:J
    for jj=1:J
        fi0(ii,jj)=(1/sqrt((pi*sigma^2)))*...
            exp(1i*x(ii)*kx+1i*y(jj)*ky...
            -(x(ii)-x0).^2/(2*sigma^2)...
            -(y(jj)-y0).^2/(2*sigma^2));
        if ((x(ii)-xv)^2+(y(jj)-yv)^2) <= a^2
            V(ii,jj)=Vmax;
        else
            V(ii,jj)=0;
        end
    end
end
n=1;
count=0;
while n < 160
    count=count+1;
    fin=solver(fi0,V,eps,eta,J);
    fi0=fin;
    n=n+1;
    fi=fin;
    [xx,yy]=meshgrid(x,y);
    colormap pink
    mesh(xx,yy,+10*V/(Vmax))
    hold
    surf(xx,yy,abs(fi))
    hold
    view(-145,75);
    axis([-0.5 1.5 -0.5 1.5 0 10])
%    colormap winter
    shading interp
    M(count)=getframe;
    drawnow;
end
end
function fin=solver(fi0,V,eps,eta,J)
lambda=2*eps^2/eta;
%
%   part 1
%
A=(1+2*1i/lambda)*diag(ones(J-2,1))-(1i/lambda)*diag(ones(J-3,1),-1)-...
    (1i/lambda)*diag(ones(J-3,1),1);
B=zeros(J-2,1);
fis=zeros(J,J);
%
%    find fis, intermediate time step.
%
for jj=2:J-1
    for ii=2:J-1
        B(ii-1)=fi0(ii,jj)+(1i/lambda)*(fi0(ii,jj+1)-2*fi0(ii,jj)+fi0(ii,jj-1));
    end
    C=linsolve(A,B);
    CC=[0;C;0];
    fis(:,jj)=CC;
end
%
%   part2
%
fiss=exp(-1i*eta*V/2).*fis;
A=(1+2*1i/lambda)*diag(ones(J-2,1))-(1i/lambda)*diag(ones(J-3,1),-1)...
    -(1i/lambda)*diag(ones(J-3,1),1);
B=zeros(J-2,1);
fins=zeros(J,J);
%
%   find new fin+1 (wavefunction at t+lambda)
%
for ii=2:J-1
    for jj=2:J-1
        B(jj-1)=fiss(ii,jj)+(1i/lambda)*(fiss(ii+1,jj)-2*fiss(ii,jj)+fiss(ii-1,jj));
    end
    C=linsolve(A,B);
    CC=[0;C;0];
    fins(ii,:)=CC.';

end
fin=exp(-1i*eta*V/2).*fins;
end