Haven’t been posting for a while, but here’s something new… Earlier I showed how to solve the 1D Schrödinger equation numerically in different situations. Now I’m going to show how to calculate the evolution of a 2D wavepacket in a potential energy field that has been constructed to mimic the classical “two-slit experiment” which shows how the mechanics of low-mass particles like electrons can exhibit interference similar to the mechanics of classical waves (sound, light, water surface, and so on).

A 2D Schrödinger equation for a single particle in a time-independent background potential V(x,y) is

Where the particle mass has been set to 1 and the Planck’s constant to .

To solve this numerically, we need the Crank-Nicolson method, as was the case when solving the 1D problem. More specifically, the linear system to be solved is

with

where the wavefunction now has two position indices and one time index, and the potential energy has only two position indices.

To form a model of the two-slit experiment, we choose a domain 0 < x < 6; 0 < y < 6 and make a potential energy function defined by

IF (x < 2.2 OR x > 3.8 OR (x > 2.7 AND x < 3.3)) THEN IF (3.7 < y < 4) THEN V(x,y) = 30

IF (x < 0.5 OR x > 5.5 OR y < 0.5 OR y > 5.5) THEN V(x,y) = 30

Otherwise V(x,y) = 0.

which corresponds to having hard walls surrounding the domain and a barrier with two holes around the line y = 3.85

For an initial condition, we choose a Gaussian wavepacket that has a nonzero expectation value of the momentum in y-direction:

An R-Code that solves this problem for a time interval 0 < t < 1 is

library(graphics) #load the graphics library needed for plotting

lx <- 6.0 #length of the computational domain (both x and y)

lt <- 1.0 #length of the simulation time interval

nx <- 47 #number of discrete lattice points in x and y directions

nt <- 60 #number of timesteps

dx <- lx/nx #length of one discrete lattice cell, same in x and y directions

dt <- lt/nt #length of timestepV = matrix(nrow=nx,ncol=nx) #potential energies at discrete xy points

for(j in c(1:nx)) { #construct the potential field with a double-slit barrier

for(k in c(1:nx)) {

V[j,k] = 0+0i

if((j*dx<2.2)||(j*dx>3.8)||((j*dx>2.7) && (j*dx<3.3))) {

if((k*dx>3.7) && (k*dx<4.0)) {

V[j,k] = 30+0i #No significant density is going to go through these barriers

}

}

if((j*dx>5.5) || (j*dx<0.5) || (k*dx>5.5) || (k*dx<0.5)) {

V[j,k] = 30+0i

}

}

}kappa1 = (1i)*dt/(2*dx*dx) #an element needed for the matrices

kappa2 <- c(1:(nx*nx)) #another elementfor(j in c(1:nx)) {

for(k in c(1:nx)) {

kappa2[(j-1)*nx+k] <- kappa1*2*dx*dx*V[j,k]

}

}psi = c(1:(nx*nx)) #array for the wave function values

for(j in c(1:nx)) {

for(k in c(1:nx)) {

psi[(j-1)*nx+k] = exp(-2*(k*dx-1.7)*(k*dx-1.7)-2*(j*dx-3)*(j*dx-3)+(2i)*k*dx) #Gaussian initial wavefunction

if((j*dx > 5.5)||(j*dx < 0.5)||(k*dx > 5.5)||(k*dx < 0.5)) {

psi[(j-1)*nx+k] = as.complex(0)

}

}

}xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points

yaxis <- c(1:nx)*dx #the y valuesA = matrix(nrow=nx*nx,ncol=nx*nx) #matrix for forward time evolution

B = matrix(nrow=nx*nx,ncol=nx*nx) #matrix for backward time evolution

P = matrix(nrow=nx,ncol=nx) #matrix for the solution after time stepping

IP = matrix(nrow=4*nx,ncol=4*nx) #matrix for the higher resolution image obtained by interpolationfor(j in c(1:(nx*nx))) { #Set the values for time evolution matrix elements

for(k in c(1:(nx*nx))) {

A[j,k]=0

B[j,k]=0

if(j==k) {

A[j,k] = 1 + 4*kappa1 + kappa2[j]

B[j,k] = 1 – 4*kappa1 – kappa2[j]

}

if((k==j+1) || (k==j-1)) {

A[j,k] = -kappa1

B[j,k] = kappa1

}

if((k==j+nx)||(k==j-nx)) {

A[j,k] = -kappa1

B[j,k] = kappa1

}

}

}for(k in c(1:nt)) { #main time stepping loop

sol <- solve(A,B%*%psi) #solve the system of equations

for(l in c(1:(nx*nx))) {

psi[l] <- sol[l]

}for(l in c(1:nx)) {

for(m in c(1:nx)) {

P[l,m] = abs(psi[(l-1)*nx + m])*abs(psi[(l-1)*nx + m]) #square of the absolute value of wave function

if(abs(V[l,m]) > 5) P[l,m] = 2

}

}

for(l in c(1:(nx-1))) {

for(m in c(1:(nx-1))) { #make a bitmap with 4 times more pixels, using linear interpolation

IP[4*l-3,4*m-3] = P[l,m]

IP[4*l-2,4*m-3] = P[l,m]+0.25*(P[l+1,m]-P[l,m])

IP[4*l-1,4*m-3] = P[l,m]+0.5*(P[l+1,m]-P[l,m])

IP[4*l,4*m-3] = P[l,m]+0.75*(P[l+1,m]-P[l,m])

}

}for(l in c(1:(4*nx))) {

for(m in c(1:(nx-1))) {

IP[l,4*m-2] = IP[l,4*m-3]+0.25*(IP[l,4*m+1]-IP[l,4*m-3])

IP[l,4*m-1] = IP[l,4*m-3]+0.5*(IP[l,4*m+1]-IP[l,4*m-3])

IP[l,4*m] = IP[l,4*m-3]+0.75*(IP[l,4*m+1]-IP[l,4*m-3])

}

}jpeg(file = paste(“plot_abs_”,k,”.jpg”,sep=””)) #save the image

image(IP, zlim = c(0,0.15))

dev.off()}

The code produces a sequence of image files, where the probability density is plotted with colors, as an output. Some representative images from this sequence (converted to grayscale) is shown below:

A video of the time evolution is shown below:

The treshold for maximum white color has been chosen to be quite low, to make the small amount of probability density that crosses the barrier visible.

The discrete grid of points has been made quite coarse here to keep the computation time reasonable, and the resolution has been increased artificially by using linear interpolation between the discrete points.

So, now we’ve seen how to solve the motion of 2D wavepackets moving around obstacles. In the next numerical methods post, I’ll go through the numerical solution of a nonlinear PDE.