Categories
Computational Physics Quantum Mechanics

Numerical Quantum Mechanics: The Time-dependent Schrödinger Equation (II)

In my previous blog post, I had an introductory look at the problem of integrating the time-dependent Schrödinger equation (TDSE) numerically, plotting the wave function at each time step, and creating an animation of the result. In this post I will describe how to solve the TDSE in more detail.

Recall the TDSE for a spinless particle in a potential V:

i\hbar \frac{\partial}{\partial t} \psi = \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \psi.

We need to do two things;

  1. Evaluate the Hamiltonian acting on the wave function, \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \psi.
  2. Integrate the TDSE itself.

Evaluating H\psi

The Hamiltonian consists of two parts; the kinetic energy operator (which is proportional to the Laplacian) and the potential operator. Both of these act on the wave function, producing new fields of complex numbers (-\frac{\hbar}{2m} {\nabla}^2 \psi and V \psi). The potential term is easy; at every point on our grid we simply multiply the wave function with the potential. The kinetic term requires us to evaluate the action of the Laplacian on the wave function, and we need to come up with a numerical scheme for this.

Spatial derivatives

The Laplacian is the sum of the (unmixed) second partial derivatives. Since we are considering a 2D space, we have:

\nabla^2 = \frac{\partial^2}{\partial x^2} +\frac{\partial^2}{\partial y^2}.

In order to evaluate the kinetic term, therefore, we must evaluate the second derivative with respect to x and y. For a good introduction to spatial derivatives in finite-difference schemes, see this Wikipedia page. A popular choice is the central difference scheme:

f'(x) \approx \frac{f(x+h) - f(x-h)}{2h},

where h is the grid spacing. A second order version also exists:

f''(x) \approx \frac{f(x+h) - 2 f(x) + f(x-h)}{h^2},

but I get better (more stable) results by simply applying the first order method twice. Assuming we pick the second order method, we end up with the following expression:

{\nabla}^2 \psi(x,y) \approx\frac{\psi(x+h, y) - 2\psi(x, y) +\psi(x-h, y)}{h^2} +\frac{\psi(x, y+h) - 2\psi(x, y) +\psi(x, y-h)}{h^2}.

Boundary conditions

Boundary conditions probably deserve a few dedicated blog posts; I will not go into them very deeply here, because in all of the videos I am presenting here, the wave function remains more or less localized in the middle of the screen, barely touching the edges. Suffice it to say for now that there are different options; one can have reflecting boundaries (so that the entire wave function remains in our simulation volume) or absorbing boundaries (so that (parts of) the wave function can propagate out of our simulation volume). Note that in the second case, the wave functions norm will generally not be preserved!

In the code I am describing in these blog posts, I use a crude scheme where I simply set \psi to zero at the boundary after each step, and only evaluate the Laplacian in the innermost region of the simulation volume (excluding a boundary zone one with a thickness of one grid point). This is not a proper implementation of boundary conditions, and I will have to revisit this issue at a later date. (Update: I’ve written a blog post about proper boundary conditions for the TDSE.)

Integrating the TDSE

In my previous post I mentioned the simplest possible explicit integration scheme for the TDSE: the forward Euler method. It works as follows:

\psi_{n+1} = \psi_{n} + \frac{\Delta t}{i \hbar} \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \psi_{n}\\ t_{n+1}=t_n+\Delta t.

While this method is easy to implement, it is a first order scheme, and the error per step is proportional to the square of the step size. This means that unless the step size is extremely small (in which case we need many integration steps, which is computationally expensive), the solution quickly becomes unstable and ‘blows up’. We therefore use the 4th order Runge-Kutta method (abbreviated RK4), which is much more stable. For more detailed information about integration schemes, I recommend the excellent Wikipedia articles on the Euler method and Runge-Kutta methods.

Let’s write out the RK4 integration scheme for the TDSE explicitly. In this scheme, the new wave function \psi_{n+1} depends not only on the current wave function \psi_{n}, but also on four ‘intermediate’ wave functions which must be computed in order:

\psi_{k_1}= \frac{1}{i \hbar} \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \psi_{n}\\ \psi_{k_2}= \frac{1}{i \hbar} \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \left( \psi_n + \frac{\Delta t}{2}  \psi_{k_1} \right) \\ \psi_{k_3}= \frac{1}{i \hbar} \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \left( \psi_n + \frac{\Delta t}{2} \psi_{k_2} \right) \\ \psi_{k_4}= \frac{1}{i \hbar} \left[-\frac{\hbar}{2m} {\nabla}^2 + V \right] \left( \psi_n + \Delta t \psi_{k_3} \right).

Having obtained the intermediate wave functions, we then compute the new wave function \psi_{n+1} as follows:

\psi_{n+1}=\psi_{n} + \frac{\Delta t}{6} \left( \psi_{k_1} + 2\psi_{k_2} + 2\psi_{k_3} + \psi_{k_4}\right)\\ t_{n+1}=t_n+\Delta t.

Wrap-up

That’s all I wanted to say about evaluating the action of the Hamiltonian on the wave function. In the next blog post, we’ll focus on boundary conditions, as well as implementing the equations discussed so far in Python, and look at some examples and plots to make things more clear. If you have any questions, don’t hesitate to post a comment and ask them!

3 replies on “Numerical Quantum Mechanics: The Time-dependent Schrödinger Equation (II)”

Hi David, looks like you are completely correct – thanks so much for pointing this out! I had it right in my code, but definitely posted the wrong formula on my blog.

Liked by 1 person

Leave a comment