Skip to main content

### Crank–Nicolson method

The Crank–Nicolson method is a popular choice when one is faced with the problem of solving the time-dependent Schrödinger equation numerically. Here, we outline the main idea of the method and apply it to a simple scenario of a one-dimensional quantum particle moving inside a box.

### Schrödinger equation

For all practical purposes, a state $$\Psi(t)$$ of a quantum system can be viewed as a complex, time-dependent vector: $\Psi(t) = \begin{pmatrix} \psi_1(t) \\ \psi_2(t) \\ \cdots \\ \psi_N(t) \end{pmatrix} \,.$ Here, $$\psi_1(t)$$, $$\psi_2(t)$$, $$\ldots$$, $$\psi_N(t)$$ are some complex-valued functions of time $$t$$. The time-evolution of $$\Psi(t)$$ is governed by the Schrödinger equation, $i \hbar \frac{d \Psi(t)}{d t} = H \Psi(t) \,,$ where $$i$$ is the imaginary unit, $$\hbar$$ is the Planck constant, and $H = \begin{pmatrix} H_{11} & H_{12} & \cdots & H_{1N} \\ H_{12}^* & H_{22} & \ldots & H_{2N} \\ \cdots & \cdots & \cdots & \cdots \\ H_{1N}^* & H_{2N}^* & \cdots & H_{NN} \\ \end{pmatrix}$ is the Hamiltonian. (Note that for $$H$$ to be a Hermitian matrix, its diagonal elements must be real.) Hereinafter, we assume that $$H$$ does not depend on time.

### Time-evolution operator

It is easy to check that the time-dependent Schrödinger equation admits the following solution: $\Psi(t) = U(t) \Psi(0) \,,$ where $U(t) = \exp \left( -\frac{i t}{\hbar} H \right) \,.$ The matrix $$U(t)$$ is called the time-evolution operator, as it evolves an initial state $$\Psi(0)$$ through time $$t$$ into the final state $$\Psi(t)$$.

The formal expression defining the evolution operator has to be understood in terms of the Taylor expansion of the exponential function, $e^z = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \ldots$ Taking $$z = (-i t / \hbar) H$$, one obtains the following infinite series representation of the evolution operator: $U(t) = I - \frac{i t}{\hbar} H - \frac{t^2}{2 \hbar^2} H^2 + \frac{i t^3}{6 \hbar^3} H^3 + \ldots$ Here $$I$$ is the $$N \times N$$ unit matrix.

### Discretization of time

In order to evolve the initial state, $$\Psi(0)$$, through time $$t$$, one discretizes $$t$$ into $$n$$ sufficiently short intervals $$\tau$$, so that $t = n \tau \,,$ and applies the evolution operator successively: \begin{align*} &\Psi(\tau) = U(\tau) \Psi(0) \\ &\Psi(2 \tau) = U(\tau) \Psi(\tau) \\ &\Psi(3 \tau) = U(\tau) \Psi(2 \tau) \\ &\cdots\cdots\cdots\cdots\cdots\cdots \\ &\Psi(t) = U(\tau) \Psi\big( (n-1) \tau \big) \end{align*} The problem is now reduced to finding an accurate and efficient numerical approximation for $$U(\tau)$$ where $$\tau$$ is small.

### Naive approach

A natural line of thinking would be to keep only the first few terms of the power series expansion defining $$U(\tau)$$. Thus, keeping only the first two terms of the expansion, one would approximate the evolution operator by $U(\tau) \simeq I - \frac{i \tau}{\hbar} H \,.$ The main drawback of this approximation is that it is not unitary: if one applied the (approximate) evolution operator twice, first forward and then backward in time, one would not recover the initial sate. Indeed, \begin{align*} U(-\tau) U(\tau) \Psi(0) &\simeq \left( I + \frac{i \tau}{\hbar} H \right) \left( I - \frac{i \tau}{\hbar} H \right) \Psi(0) \\ &= \left( I + \frac{\tau^2}{\hbar^2} H^2 \right) \Psi(0) \\ &\not= \Psi(0) \,. \end{align*} At each evolution step, this approximation introduces an error to the norm of the state. For instance, the norm of the (approximated) state $$U(\tau) \Psi(0)$$ is not the same as the norm of $$\Psi(0)$$. The error in the norm grows over time and renders this naive method unstable. Increasing the number of terms retained in the power series expansion of $$U(\tau)$$ does not resolve the issue.

### Crank–Nicolson approach

The Crank–Nicolson approximation of $$U(\tau)$$ is based on approximating the exponential function by a simple rational function (the so-called $$[1/1]$$ Padé approximation), $e^z \simeq \frac{1 + z/2}{1 - z/2} \,,$ valid for sufficiently small $$z$$. The validity of this approximation is evident from the Taylor expansion of the right-hand side around zero: $\frac{1 + z/2}{1 - z/2} = 1 + z + \frac{z^2}{2} + O(z^3) \,.$ Here $$O(z^3)$$ represents terms of order $$z^3$$ and higher.

Using $$z = (-i\tau/\hbar) H$$, one obtains the Crank-Nicolson approximation of the evolution operator: $U(\tau) \simeq \left(I + \frac{i \tau}{2 \hbar} H \right)^{-1} \left(I - \frac{i \tau}{2 \hbar} H \right) \,.$ The main advantage of this approximation, compared to the naive approach discussed above, is that it respects the unitarity of the evolution operator, $U(-\tau) U(\tau) = I \,,$ and hence preserves the norm of the quantum state.

### Time-evolution step

According to the Crank–Nicolson approximation, the relation between an (unknown) state of the system at time $$(k+1) \tau$$ and a (known) state at time $$k \tau$$ is \begin{align*} \Psi\big( (k+1) \tau \big) &= U(\tau) \Psi(k \tau) \\ &\simeq \left(I + \frac{i \tau}{2 \hbar} H \right)^{-1} \left(I - \frac{i \tau}{2 \hbar} H \right) \Psi(k \tau) \,. \end{align*} The relation can also be written as $\left(I + \frac{i \tau}{2 \hbar} H \right) \Psi\big( (k+1) \tau \big) \simeq \left(I - \frac{i \tau}{2 \hbar} H \right) \Psi(k \tau) \,.$ This means that the sought state $$\Psi\big( (k+1) \tau \big)$$ is the solution vector $$\boldsymbol{x}$$ of following system of linear equations: $A \boldsymbol{x} = \boldsymbol{b} \,,$ where $$A$$ is a matrix defined as $A = I + \frac{i \tau}{2 \hbar} H$ and $$\boldsymbol{b}$$ is a vector given by $\boldsymbol{b} = \left(I - \frac{i \tau}{2 \hbar} H \right) \Psi(k \tau) \,.$ This way, the problem of propagating a quantum state through $$\tau$$ has been reduced to the problem of solving a system of linear equations. The latter is standard, and appropriate library functions are readily available in any major programming language.

### An example in Python

The Python code below uses the Crank–Nicolson method to simulate the motion of a one-dimensional quantum particle between two walls at $$x=0$$ and $$x=1$$. Hereinafter we use the atomic units in which the Planck constant and particle’s mass are both equal to unity.

import numpy as np
from scipy.linalg import solve_banded
import matplotlib.pyplot as plt

Nx = 1000  # number of spatial points
x  = np.linspace(0, 1, Nx)  # spatial interval
dx = x[1]-x[0]  # spatial step

Nt = 1000  # number of time steps
dt = 0.000005  # time step

psi  = np.exp(-200*(x-0.2)**2 + 1j*100*x)  # Gaussian wave packet
norm = np.trapz(np.abs(psi)**2, x=x)
psi  = psi / np.sqrt(norm)  # normalized initial wave function

A = 0j*np.zeros((3,Nx))  # matrix (I+H*dt/2) encoded for 'solve_banded' method
for i in range(Nx):
A[0,i] = -1j * (dt/2) / (2*dx**2)
A[1,i] = 1 + 1j * (dt/2) / dx**2
A[2,i] = -1j * (dt/2) / (2*dx**2)

fig, ax = plt.subplots()
ax.axis('off')

lw = 1  # line thickness
ax.plot(x, np.abs(psi)**2, linewidth=lw)  # plot initial probability density
plt.gcf().text(0.13, 0.9, 'INITIAL', fontsize=12, color='gray')

for k in range(Nt):  # loop over time

Hpsi = 0j*np.zeros(Nx)  # Hamiltonian H=(-1/2)*d^2/dx^2 applied to psi
Hpsi[0] = (2*psi[0] - psi[1]) / (2*dx**2)
for i in range(1,Nx-2):
Hpsi[i] = (-psi[i-1] + 2*psi[i] - psi[i+1]) / (2*dx**2)
Hpsi[Nx-1] = (-psi[Nx-2] + 2*psi[Nx-1]) / (2*dx**2)

b = psi - 1j * (dt/2) * Hpsi  # matrix (I-H*dt/2) applied to psi

psi = solve_banded((1,1), A, b)  # solve (I+H*dt/2) psi_new = (I-H*dt/2) psi

if (k+1) % 100 == 0:  # plot probability density every 100 time steps
lw += 0.25
ax.plot(x, np.abs(psi)**2, linewidth=lw)

plt.gcf().text(0.72, 0.41, 'FINAL', fontsize=12, color='gray')

plt.tight_layout()
plt.show()


The code produces the following output: