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.

The source code for the double slit for schrodinger 2d barrier is not ok.

LikeLike

Where is the problem, in the Crank-Nicolson finite differencing or somewhere else? I wrote this three years ago, so I can’t remember whether I checked the conservation of norm of the wave function and so on.

LikeLike

Ok, now I noticed it after looking again, there’s something crazy happened when copying that from the source code file. I’ll fix that.

LikeLike

Thanks, but how i should simulate the repulsion or atraction from a EM field in the potential matrix.

I mean, if I have the same charge the repulsive forze is like the potential barrier that stops, but is the charge is one postive and the oher negative, it have to atract the wave function, but putting the potential energy negative it will do the same. ¿ or not?

LikeLike

If it’s a static electric field, then it can be modelled with a time-independent potential energy V(x,y) but this ignores radiation emission by the charged particle at high acceleration. If you have a positive and negative charged particle interacting in 2 dimensions, it will be a 4-dimensional problem due to both particles having an x and y coordinate. A Crank-Nicolson time stepping of a 4D wave function usually takes too much time to be done on a home computer.

LikeLike

Thanks, but if I replace your barrier with an electron (or proton)I and fill the potential matrix in your code with the values at certain distance in fm from a EM field ( also I will add strong forze too ):like this:

EMPotFromDist(double distMeters){return (2.30668E-28/distMeters)*JulestoMev;}

It will be ok? But I think this is time dependant,, because if I calc the time independant, for each step I have to move the particle in space getting the speed from the E kinetic (and it changes with the EM pot),.. Or not?

And my question before was, what else i have to change in your code to make the barrier atractive, (the potential,,the inital gaussian wavepacket)?.

LikeLike

You can’t represent an interaction with another particle as a static potential energy V(x,y), and not even a time-dependent V(x,y,t) unless the other particle is much heavier and follows some known trajectory.

Making the barrier attractive requires nothing more than changing the sign of function V(x,y).

LikeLike

Thanks, but I changed the sign of V(x,y) and it still repulsive. .in

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

I have another example of Hydrogen atom, and it is time independant using n,l,m,but i would like to see how the orbit change in time with a gaussian wave packet like yours. With yours I put 2 repulsive particles with a potential spherical probability using the EM forze by distance, placed on the mass center of the probability, but to simulate atraction i need to know how to model this potential or wave equation

Also I translate to phyton embended in c++and use sparse matrix with spsolve(0.02 secs) instead of solve (8 secs) because is much faster.

LikeLike

ok, thanks in this link there is source code for 2 particle scatering guassian animation. They calc in 2d when its a 1d correlationm, and then extract the 1d solution for both

LikeLike