Last Updated on November 30, 2025 by Dr. Sushanta Barman
Numerical simulation of wave packet dynamics plays an important role in understanding quantum systems. It helps us study how particles behave at very small scales. These simulations support the development of tools and technologies used in quantum optics, chemical physics, and nanotechnology.
Spectral methods, such as the split-step Fourier method, work very well when the Hamiltonian can be separated into kinetic and potential energy terms. However, this is not always possible for complex potentials. In such cases, finite-difference methods like the Crank–Nicolson (CN) scheme are very useful. CN is unconditionally stable and provides high accuracy for a wide range of potentials.
The Crank–Nicolson method was introduced by John Crank and Phyllis Nicolson in 1947 for solving the heat equation [1]. It extended earlier implicit finite-difference ideas by providing stability and second-order accuracy in time. Since then, it has become a widely used numerical scheme for partial differential equations in physics and engineering.
In this tutorial, we will derive the CN discretization of the one-dimensional time-dependent Schrödinger equation (TDSE), show how it leads to a tridiagonal linear system, and implement the Thomas algorithm for efficient time stepping.
Finally, we will go through an example using Python code to propagate a Gaussian wave packet towards a square potential barrier and calculate its transmission probability.
The Time-Dependent Schrödinger Equation
The time-dependent schrödinger equation (TDSE) in one dimension is given by,
\[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),\,\,\,\,\,(1)\]
where \(\psi(x,t)\) denotes the wave function, \(m\) is the particle mass, \(\hbar\) is the reduced Planck constant, and \(V(x)\) represents the potential energy.
The TDSE is also written as,
\[ i\hbar\,\frac{\partial \psi(x,t)}{\partial t} = H\psi(x,t),\,\,\,\,\,(2)\]
where \( H \) is the Hamiltonian operator given by, \( H = -\frac{\hbar^2}{2m}\,\frac{\partial^2}{\partial x^2} + V(x). \)
We seek a stable, accurate scheme to march \(\psi\) forward in time by discrete steps \(\Delta t\).
The Crank–Nicolson Method
Idea behind CN scheme
In this method, we estimate the derivative at the midpoint in time as,
\[\frac{\partial \psi}{\partial t} \approx \frac{\psi^{n+1}-\psi^n}{\Delta t}.\]
But the right-hand side uses an average of Hamiltonian terms:
\[H\psi \approx \frac{1}{2}(H\psi^{n+1} + H\psi^n)\]
Inserting it into the TDSE (Eq. (2)),
\[i\hbar \frac{\psi^{n+1}-\psi^n}{\Delta t} = \frac{1}{2}H(\psi^{n+1} + \psi^n)\]
After rearranging,
\[\underbrace{\left(I + \frac{i\Delta t}{2\hbar}H\right)}_{A} \psi^{n+1} = \underbrace{\left(I – \frac{i\Delta t}{2\hbar}H\right)}_{B} \psi^n\]
This gives a linear equation,
\[A\psi^{n+1} = B\psi^n.\,\,\,\,\,(3)\]
The matrices \(A\) and \(B\) are tridiagonal (due to second-order spatial finite-difference), so solving is very fast!
To advance the wave function from time step \(n\) to \(n+1\), the Crank–Nicolson scheme solves the above linear matrix equation (Eq. (3)).
Explicit form of the linear system
We first define the constants:
\[\alpha = \frac{i\hbar \Delta t}{4m \Delta x^2}, \qquad \beta = \frac{i \Delta t}{2\hbar}, \qquad t_0 = \frac{\hbar^2}{m\Delta x^2} \]
Matrix \(A\): Implicit operator (left-hand side) of Eq. (3)
The matrix \(A\) contains the future-time wave function values and ensures numerical stability. Its entries are
\[A_{j,j} = 1 + \beta\left(t_0 + V_j\right), \qquad A_{j,j\pm1} = -\alpha\]
Thus, \(A\) is a tridiagonal matrix,
\[A =\begin{pmatrix}1+\beta(t_0+V_1) & -\alpha & 0 & \cdots & 0 \\-\alpha & 1+\beta(t_0+V_2) & -\alpha & \cdots & 0 \\0 & -\alpha & 1+\beta(t_0+V_3) & \cdots & 0 \\\vdots & \vdots & \vdots & \ddots & -\alpha \\0 & 0 & 0 & -\alpha & 1+\beta(t_0+V_{N-1})\end{pmatrix}\]
Matrix \(B\): Explicit operator (right-hand side) of Eq. (3)
The matrix \(B\) contains values from the current time step and plays the role of a prediction step. Its entries are
\[ B_{j,j} = 1 – \beta\left(t_0 + V_j\right), \qquad B_{j,j\pm1} = +\alpha \]
The matrix form is
\[B =\begin{pmatrix}1-\beta(t_0+V_1) & +\alpha & 0 & \cdots & 0 \\+\alpha & 1-\beta(t_0+V_2) & +\alpha & \cdots & 0 \\0 & +\alpha & 1-\beta(t_0+V_3) & \cdots & 0 \\\vdots & \vdots & \vdots & \ddots & +\alpha \\0 & 0 & 0 & +\alpha & 1-\beta(t_0+V_{N-1})\end{pmatrix}\]
Finally, Eq. (3) must be solved at every time step to obtain \(\psi^{n+1}\). Because \(A\) is tridiagonal, the system can be efficiently solved using the Thomas algorithm, which scales as \(O(N)\).
Boundary Conditions
To enforce hard-wall boundaries (infinite potential well condition):
\[ \psi_0 = 0, \qquad \psi_N = 0\]
Thus, the first and last components are fixed and are not included in the solved system.
Numerical Solution via the Thomas Algorithm
Once the Crank-Nicolson discretization is written in matrix form (Eq. (3)), our task at each time step is to solve for the new wave function \(\psi^{n+1}\).
Directly inverting \(A\) would require \(\mathcal{O}(N^3)\) computational effort, which becomes extremely expensive for realistic grid sizes.
However, the tridiagonal structure allows us to use the Thomas algorithm, a simplified version of Gaussian elimination that reduces the cost to only \(\mathcal{O}(N)\).
This is what makes the CN method fast enough for real quantum simulations.
What the Thomas algorithm does
The equation,
\[A \psi^{n+1} = d, \]
where \(d = B\psi^n\) represents a system of (N-1) linear equations in (N-1) unknowns (excluding boundary nodes). The Thomas algorithm solves this in two phases:
(1) Forward sweep: Elimination
We eliminate the lower diagonal elements of \(A\), creating modified coefficients:
\[c’_j, \text{ and } d’_j\]
These represent the updated upper-diagonal and the updated right-hand side values after elimination. All operations move from left to right (from index \(j=1\) to \(N-1\)).
This step transforms the system into an upper triangular matrix, making back-substitution possible.
(2) Backward substitution
Once elimination is complete in the first step, we solve for the wave function elements,
\[\psi^{n+1}_{N-1},\ \psi^{n+1}_{N-2},\ \dots,\ \psi^{n+1}_{1}\]
All operations move from right to left. This phase efficiently reconstructs \(\psi^{n+1}\) without explicitly inverting \(A\).
Why we precompute coefficients
The matrix \(A\) does not change in time for a time-independent potential \(V(x)\).
Therefore, the forward-sweep coefficients \(c’_j\) can be computed once before the time-loop begins. This significantly speeds up the full simulation.
Full CN + Thomas update at each time step
\[\bbox[5px,border:1px solid #000]{\text{Given }\psi^n \Rightarrow \begin{cases}d = B \psi^n & \text{(construct RHS)} \\A \psi^{n+1} = d & \text{(solve using Thomas algorithm)}\end{cases}}\]
Finally, enforce Dirichlet boundary conditions,
\[\psi_0^{n+1} = \psi_N^{n+1} = 0.\]
Key advantages for CN
The Crank–Nicolson method offers several key computational advantages: (i) it has a linear computation cost of \(\mathcal{O}(N)\), (ii) it avoids explicit matrix inversion, which makes the algorithm more stable and efficient, and (iii) it scales well for large spatial grids, typically in the range of \(10^5 – 10^6\) points.
Description of the Simulated Physical System
The numerical example considers a free electron wave packet traveling toward a square potential barrier located at the center of a one-dimensional spatial domain, as shown in Fig. 1.
The initial state is a normalized Gaussian wave packet with width \(\sigma_x = 20\ \text{nm}\), centered at \(x_0 = -200\ \text{nm}\), and assigned a positive momentum \(p_0 = \hbar k_0\), corresponding to a de Broglie wavelength of \(5\ \text{nm}\). This ensures the packet moves from left to right toward the barrier.
The potential consists of a finite-width rectangular barrier of width \(10\ \text{nm}\) and height \(1.02,(\hbar^2 k_0^2 / 2m)\), as shown in Fig. 1. Since the barrier height is slightly above the kinetic energy of the electron, the wave packet undergoes quantum tunneling, resulting in partial transmission and partial reflection.
Hard-wall boundary conditions \(\psi = 0\) at both edges mimic an infinite potential well that prevents the wave packet from leaving the simulation domain.
Python Code
Below is a complete Python implementation. Comments highlight key steps.
📄 Show Code
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import FFMpegWriter
# --------------------------------------------
# Global style for plots
# --------------------------------------------
plt.rcParams.update({
"font.family": "serif",
"mathtext.fontset": "stix",
"font.size": 16,
"axes.labelsize": 16,
"axes.titlesize": 16,
"axes.linewidth": 1.2,
"xtick.direction": "in",
"ytick.direction": "in",
"xtick.top": True,
"ytick.right": True,
"xtick.major.size": 6,
"ytick.major.size": 6,
"legend.frameon": False
})
# --------------------------------------------
# 1. Physical constants and basic parameters
# --------------------------------------------
hbar = 1.0545718e-34
m_e = 9.10938356e-31
wavelength = 5e-9
k0 = 2 * np.pi / wavelength
# --------------------------------------------
# 2. Spatial and temporal grids
# --------------------------------------------
Lx = 1e-6
Nx = 4024
dx = Lx / Nx
x = np.linspace(-Lx/2, Lx/2, Nx)
dt = 0.2e-15
Nt = 16000
# --------------------------------------------
# 3. Initial Gaussian wave packet
# --------------------------------------------
sigma_x = 20e-9
x0 = -2e-7
kx0 = k0
def initial_wavefunction():
psi_init = np.exp(-(x - x0)**2 / (2 * sigma_x**2)) * np.exp(1j * kx0 * x)
norm_init = np.sqrt(np.sum(np.abs(psi_init)**2) * dx)
psi_init /= norm_init
return psi_init
psi = initial_wavefunction()
# --------------------------------------------
# 4. Potential barrier
# --------------------------------------------
barrier_width = 10e-9
barrier_height = 1.02 * hbar**2 * k0**2 / (2 * m_e)
V = np.zeros_like(x)
V[np.abs(x) < barrier_width / 2] = barrier_height
# --------------------------------------------
# 5. Crank–Nicolson
# --------------------------------------------
alpha = 1j * hbar * dt / (4 * m_e * dx**2)
beta = 1j * dt / (2 * hbar)
t0 = hbar**2 / (m_e * dx**2)
off_diag_A = -alpha
off_diag_B = alpha
diag_A = 1 + beta * (t0 + V)
diag_B = 1 - beta * (t0 + V)
a = np.full(Nx - 1, off_diag_A, dtype=complex)
b = diag_A.copy()
c = np.full(Nx - 1, off_diag_A, dtype=complex)
c_prime = np.zeros(Nx - 1, dtype=complex)
c_prime[0] = c[0] / b[0]
for i in range(1, Nx - 1):
denom = b[i] - a[i - 1] * c_prime[i - 1]
c_prime[i] = c[i] / denom
def thomas_solve(d):
d_prime = np.zeros_like(d, dtype=complex)
d_prime[0] = d[0] / b[0]
for i in range(1, Nx):
denom = b[i] - a[i - 1] * c_prime[i - 1]
d_prime[i] = (d[i] - a[i - 1] * d_prime[i - 1]) / denom
psi_new = np.zeros_like(d, dtype=complex)
psi_new[-1] = d_prime[-1]
for i in range(Nx - 2, -1, -1):
psi_new[i] = d_prime[i] - c_prime[i] * psi_new[i + 1]
return psi_new
# --------------------------------------------
# 7. Animation setup
# --------------------------------------------
fig, ax = plt.subplots(figsize=(8, 4.5))
rho = np.abs(psi)**2
rho_line, = ax.plot(
x * 1e9,
rho / np.max(rho),
color='b',
lw=2.0,
label=r"$|\psi(x,t)|^2$"
)
V_scaled = V / np.max(V) if np.max(V) != 0 else V
V_line, = ax.plot(
x * 1e9,
V_scaled,
color='k',
linestyle="-",
lw=2.0,
label=r"$V(x)$ (scaled)"
)
ax.set_xlabel(r"$x\ \mathrm{(nm)}$")
ax.set_ylabel(r"Probability density (arb. units)")
title = ax.set_title(f"Step 0 / {Nt} Transmission = 0.0000")
ax.legend(loc="upper right")
ax.set_ylim(-0.05, 1.5)
fig.tight_layout()
fig.savefig("wavepacket_initial_t0.png", dpi=400, bbox_inches="tight")
print("Initial plot saved as: wavepacket_initial_t0.png")
steps_per_frame = 50
n_blocks = Nt // steps_per_frame
n_frames = n_blocks + 1
def update(frame):
global psi
if frame == 0:
psi = initial_wavefunction()
current_step = 0
else:
for _ in range(steps_per_frame):
d = diag_B * psi
d[1:] += off_diag_B * psi[:-1]
d[:-1] += off_diag_B * psi[1:]
psi = thomas_solve(d)
psi[0] = 0.0 + 0.0j
psi[-1] = 0.0 + 0.0j
norm = np.sqrt(np.sum(np.abs(psi)**2) * dx)
psi /= norm
current_step = frame * steps_per_frame
if current_step > Nt:
current_step = Nt
rho = np.abs(psi)**2
transmission = np.sum(rho[x > barrier_width / 2]) * dx
rho_line.set_ydata(rho / np.max(rho))
title.set_text(f"Step {current_step} / {Nt} Transmission = {transmission:.4f}")
return rho_line, V_line, title
anim = FuncAnimation(
fig,
update,
frames=n_frames,
interval=100,
blit=False,
repeat=False
)
fig.tight_layout()
plt.show()
# --------------------------------------------
# 8. Final figure
# --------------------------------------------
rho_final = np.abs(psi)**2
fig_final, ax_final = plt.subplots(figsize=(4, 2.25))
ax_final.plot(
x * 1e9,
rho_final / np.max(rho_final),
lw=2.0,
color='b',
label=r"$|\psi(x,t_{\mathrm{final}})|^2$"
)
ax_final.plot(
x * 1e9,
V_scaled,
color='k',
linestyle="-",
lw=2,
label=r"$V(x)$ (scaled)"
)
ax_final.set_xlabel(r"$x\ \mathrm{(nm)}$")
ax_final.set_ylabel(r"Probability density (arb. units)")
ax_final.set_title(rf"Final state after {Nt} time steps")
ax_final.legend(loc="upper right")
ax_final.set_ylim(-0.05, 1.5)
fig_final.tight_layout()
fig_final.savefig(f"wavepacket_final_t{Nt}.png", dpi=400, bbox_inches="tight")
print(f"Final plot saved as: wavepacket_final_t{Nt}.png")
# --------------------------------------------
# 9. Save animation
# --------------------------------------------
writer = FFMpegWriter(fps=30)
anim.save("cn_wavepacket_tunneling.mp4", writer=writer)
print("Animation saved as: cn_wavepacket_tunneling.mp4")
| Parameter | Value |
|---|---|
| Domain length (\(L_x\)) | 1 μm |
| Number of grid points (\(N_x\)) | 4024 |
| Spatial resolution (\(\Delta x\)) | ≈ 0.25 nm |
| Time step (\(\Delta t\)) | 0.2 fs |
| Total number of steps (\(N_t\)) | 16000 |
| Initial position (\(x_0\)) | −200 nm |
| Wave-packet width (\(\sigma_x\)) | 20 nm |
| Wave number (\(k_0\)) | \(2\pi/5\,\mathrm{nm}^{-1}\) |
| Barrier width | 10 nm |
| Barrier height | \(1.02\,(\hbar^2 k_0^2 / 2m)\) |
Results and Discussion
The simulation result is presented as an animation in Video 1. The Crank–Nicolson method smoothly evolves the electron wave packet as it moves toward the potential barrier. The results show two clear quantum behaviors: wave-packet spreading and tunneling.
At the start, the Gaussian wave packet travels freely from the left. It moves toward the barrier due to its initial momentum. When the wave packet reaches the barrier, the interaction divides it into two parts:
(a) Reflection
A portion of the wave packet is reflected back to the left, as shown in Video 1. The reflected part spreads more due to the interference between the incoming and reflected waves.
(b) Transmission (Quantum Tunneling)
Another part of the wave packet tunnels through the barrier and appears on the right side, as shown in Video 1. Its amplitude is smaller because the wave decays while crossing the barrier. This confirms that the particle has less energy than the barrier height, but still has a finite probability of crossing it.
During the simulation, the total probability remains constant, showing that the Crank–Nicolson method preserves normalization. The wave-function density evolves smoothly without numerical instability or unwanted oscillations. This demonstrates that the method is stable, accurate, and suitable for simulating quantum tunneling problems.
Conclusion
The Crank–Nicolson scheme provides a robust, accurate, and stable framework for simulating time-dependent quantum dynamics with arbitrary potentials. By reducing each time‑step update to solving a simple tridiagonal linear system, one achieves \(O(N)\) performance and unconditional stability.
This tutorial walked through the mathematical derivation, implementation via the Thomas algorithm, and a complete Python code example that tracks a Gaussian wave packet interacting with a square barrier.
You can modify the code to try different potentials or simulate higher-dimensional systems to study various quantum effects.
References
- J. Crank, P. Nicolson, “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” Mathematical Proceedings of the Cambridge Philosophical Society 43, 50-67 (1947). doi:10.1017/S0305004100023197
- Sushanta Barman, “Why Physicists Love Gaussian Wave Packets” MatterWaveX.Com, May 25, (2025).
- Sushanta Barman, “Split-Step Fourier Method for TDSE: Step-by-Step Guide” MatterWaveX.Com, August 18, (2025).
- Sushanta Barman, “Schrödinger Equation for Matter-Wave dynamics” MatterWaveX.Com, September 7, (2024).
I am a Fellow of Academic and Research Excellence (FARE) in the Department of Physics at IIT Kanpur. My work focuses on ion-beam optics and matter-wave phenomena. I am also interested in emerging matter-wave technologies.
