## Double Slit Bohmian Trajectories

One of the archetypal examples of quantum mechanics is the famous double slit experiment. In this experiment, particles are individually shot towards a double slit, and even with individual massive particles, an interference pattern emerges, indicating the "particle-wave duality" of matter. More details on Wikipedia, as usual.

One can calculate the interference patterns by hand using fairly simple geometry. In this article, though, I'll show how to do it by solving the time-dependent Schrödinger equation on a computer, and how to obtain Bohmian particle trajectories.

### Schrödinger equation

Let us recall to mind the basic facts. The time-dependent Schrödinger equation is \begin{align} i\hbar \frac{\partial}{\partial t}\psi (\vec{x},t) = \bigg[\frac{-\hbar ^2}{2m}\nabla ^2 + V(\vec{x})\bigg]\psi (\vec{x},t) \end{align}

In other words, there are derivatives in two variables: position and time. Unfortunately, $$\psi$$ is a continuous function, and computers can't really represent continuous functions; we must define it on a grid. We suppose that our (for this purpose 2-dimensional) space consists of a limited number of squares of area $$(dx)^2$$. Then we only need to computer our wave function on a limited number of grid points, a task for which computers are ideally suited.

Let's first consider the right hand side of the equation. What is $$V(\vec{x})$$? Well, that must be the potential of the double slit wall that we want to model. This wall is supposed to be impenetrable, or at least close to impenetrable. High potential means that particles are unlikely to be found in the area. Therefore, we should make our potental a very large number where the wall is supposed to be, and zero everywhere else, something like this: This part is easy to compute, since you only need to do it once - just create a 2D array and set the parts corresponding to the walls to a large number.

What about the derivative? Well, the derivative denotes the change of $$\psi$$ over space, and is defined as

\begin{align} \frac{\partial \psi }{\partial x} = \lim_{h\rightarrow 0}\frac{\psi(x+h)-\psi(x)}{h} \end{align}

Note that in the general definition we divide the difference of $$\psi$$ between two points separated by $$h$$ with the distance between the points and that distance approaches the limit of 0. But, in our discrete grid, there is really no such thing as a "zero distance" since it consists of a finite number of squares. Instead, we define \begin{align} \frac{\partial \psi}{\partial x} = \frac{\psi _{n+1}-\psi_{n}}{dx} \end{align}

where $$\psi _n$$ means value of $$\psi$$ at the nth grid square. That is, we take the difference between neighbouring squares and divide by the length of the side of the squares, which is taken to be the distance between two adjacent grid points.

We need the Laplacian, $$\nabla ^2 = \frac{\partial ^2}{\partial x^2} + \frac{\partial ^2}{\partial y^2}$$. For the second derivative, we get (by applying the derivative to the first derivative again!)

\begin{align} \nabla ^2\psi_{i,j} = \frac{\psi_{i+1,j}+\psi_{i-1,j}+\psi_{i,j+1}+\psi_{i,j-1}-4\psi_{i,j}}{dx^2} \end{align}

So, as an end result, we get a bunch of equations (one for each grid point) which define the right side of the equation - and on the left side, we have the time evolution.

To evolve the system in time, we apply a numerical integrator, such as the Runge-Kutta method. More about that in the when I show some code.

### Bohmian trajectories

Now, given the wave function at time t, it is possible to calculate the Bohmian trajectories from de Broglie's guidance equation,

\begin{align} \frac{dQ(t)}{dt} = \frac{\hbar}{m}\mathrm{Im}{\frac{\nabla \psi}{\psi}}(Q,t) \end{align}

That is, we calculate the gradient of the wave function at the particle's position Q, divide by the wave function at that point, take the imaginary part of this and then advance the solution with a numerical integrator. Note that this is basically independent from the Schrödinger equation: to get the correct predictions for the double slit experiment, you need not assume the existence of Bohmian trajectories at all - but they're fun to look at regardless.

### Implementation

If you have numpy installed, you can follow along with this, but I also have a code here which contains an optimized version with all the details.

Firstly, the Laplacian can be computed by using something like this, assuming $$psi$$ is a 2d array:


def laplace(psi,dx):
deriv = np.zeros((len(psi[:,0]),len(psi[:,0])))
for i in range(1,len(psi[:,0])-1):
for j in range(1,len(psi[:,0])-1):
deriv[i,j] = (psi[i+1,j]+psi[i-1,j]+psi[i,j-1]+psi[i,j+1]-4*psi[i,j])/dx/dx
return deriv


The Runge-Kutta steps can be taken with this (dt size of the timestep, and you must precompute V):


def take_step(psi,dt)
k1 = -(1. / (2.j)) * laplace(psi,Npoints,dx) + (1. / 1.j)*V*(psi)
k2 = -(1. / (2.j)) * laplace(psi + 0.5 * dt * k1,dx)+ (1. / 1.j)*V*(psi + 0.5 * dt * k1)
k3 = -(1. / (2.j)) * laplace(psi + 0.5 * dt * k2,dx)+ (1. / 1.j)*V*(psi + 0.5 * dt * k2)
k4 = -(1. / (2.j)) * laplace(psi + dt * k3,dx)+ (1. / 1.j)*V*(psi + dt * k3)

psinew = psi + dt * 1. / 6. * (k1 + 2. * k2 + 2. * k3 + k4)
return psinew



Note that in Python, this means basically going through all the grid points at once, so k1, k2, k3 and k4 are 2D arrays of the same size as $$psi$$.

There is one more detail that is important to point out: what is $$\psi$$ at the beginning? A good choice is to use a Gaussian, as they nicely satisfy e.g. the uncertainty relation. To construct an initial Gaussian, we can use:


x = np.linspace(-1,1,201)
y = np.linspace(-1,1,201)
kx = -1000.0
sigma = 1/6 #You can test different values here.
for i in range(Npoints):
for j in range(Npoints):
psi[i,j] = np.exp(-((x[i]+0.2)**2+(y[j])**2) / 2 / sigma / sigma)*(np.cos(kx*x[i])+1j*np.sin(kx*x[i]))


This would create a Gaussian at point (-0.2,0) in our grid. More details in the github code.

The end result, when animated, looks like this (for two initial particle positions)

That's it for now. Here's the main site, and a twitter.