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: