Introduction
In this post we simulate the schrodinger equation using FDTD (finite-difference time-domain). We scatter a gaussian wavepacket off a step potential.
This so far has been one of my favorite simulations. Normally we can only exactly solve the Schrodinger equations for a few simple potentials, but this method allows us to see what would happen in any potential!
Results
Video:
Details
This work was heavily based off these two resources:
- https://www.youtube.com/watch?v=S-6Z8N-30AU
- https://scipy-cookbook.readthedocs.io/_downloads/bca7b26f2efcc85a15a70187c3989aa1/Schrodinger_FDTD.pdf
The idea is relatively simple. We want to simulate the Schrodinger equation for a specific potential:
$$ i \hbar \frac{\partial \psi(x, t)}{\partial t} = - \frac{\hbar^2}{2m} \frac{\partial^2 \psi(x, t)}{\partial x^2} + V(x) \psi(x,t)$$
Specifically we want to see how \( \psi \) changes with time.
We see that the time derivative can be calculated by the Laplacian (double derivative with respect to position) and the potential.
The trick we use is to turn the continuous equation into a lattice of points, and use the definition of a derivative:
$$ \frac{d \psi}{d x} \approx \lim_{a \to small} \frac{ \psi(x+a) - \psi(x) }{a} $$
$$ \frac{d^2 \psi}{d x^2} \approx \lim_{a \to small} \frac{ \psi(x+a) - 2\psi(x) + \psi(x-a) }{a^2} $$
Plugging this into the schrodinger equation and massaging the constants we get these two very simple update equations (for each grid spot)
this.psi_i_future[x] = this.c1 * (this.psi_r_present[x + 1] - 2 * this.psi_r_present[x] + this.psi_r_present[x - 1]) - this.c2V[x] * this.psi_r_present[x] + this.psi_i_past[x];
this.psi_r_future[x] = -this.c1 * (this.psi_i_present[x + 1] - 2 * this.psi_i_present[x] + this.psi_i_present[x - 1]) + this.c2V[x] * this.psi_i_present[x] + this.psi_r_past[x];
And that’s it!
I just love how much easier it is to code solutions numerically than attempt to solve them exactly.
Random Notes
- The article recommended splitting the schrodinger equation into two equations for the real and complex values. I followed the advice, but I think for 2D I’m just going to use complex numbers instead
- The simulation was very slow at first, thankfully Chrome has really nifty profiling tools and I was able to see I just had a heap leak and garbage collection was being called a lot
- I am slightly nervous about the 2D simulation being slow, but we’ll cross that bridge later
- The initial gaussian pulse was initially being split into two parts going opposite directions. There was a bug in my update equation which caused that to happen
Future
Up next is scattering electrons off a double slit. So I’m working on a 2D simulation of this.