Crank–Nicolson Method for the Time-Dependent Schrödinger Equation: Step-by-Step Guide

By /

Last Updated on July 4, 2025 by Sushanta Barman

Numerical simulation of quantum wave packet dynamics is at the heart of many applications in quantum optics, chemical physics, and nanotechnology. While spectral (split-step) methods shine when the Hamiltonian splits neatly into kinetic and potential parts, fully implicit finite-difference schemes like Crank–Nicolson (CN) offer unconditional stability and excellent accuracy for arbitrary potentials.

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 walk through a complete Python example: propagating a Gaussian wave packet against a square potential barrier and computing its transmission probability.

The Time-Dependent Schrödinger Equation

The 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)\,,\]

where \(\psi(x,t)\) denotes the wave function, \(m\) the particle mass, \(\hbar\) the reduced Planck constant, and \(V(x)\) the external potential.

We seek a stable, accurate scheme to march \(\psi\) forward in time by discrete steps \(\Delta t\).

Crank–Nicolson Discretization

The Crank–Nicolson method is obtained by applying the trapezoidal (midpoint) rule in time and central differences in space:

  1. Time-center the derivative: \[\frac{\psi^{n+1}_j – \psi^n_j}{\Delta t}\;\approx\;\left.\frac{\partial\psi}{\partial t}\right|_{t = t_n + \frac{\Delta t}{2}}.\]
  2. Average the right-hand side at times \(t_n\) and \(t_{n+1}\): \[i\hbar\,\frac{\psi^{n+1}_j – \psi^n_j}{\Delta t}=\frac{1}{2}\Bigl[\hat H\,\psi^{n+1}_j + \hat H\,\psi^n_j\Bigr],\] where \(\hat H = -\frac{\hbar^2}{2m}\partial_{xx} + V_j\).
  3. Discretize second derivative by central differences: \[\partial_{xx}\psi_j \approx \frac{\psi_{j+1} – 2\psi_j + \psi_{j-1}}{\Delta x^2}.\]

Putting it all together leads to the linear system,
\[A\,\psi^{\,n+1} = B\,\psi^n,\]
with tri-diagonal matrices \(A\) and \(B\).

Form of the Matrices

Define

\[\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}.\]

Then for grid point \(j\):

  • Off-diagonal entries of both \(A\) and \(B\): \[(A)_{j,j\pm1} = -\,\alpha,\quad (B)_{j,j\pm1} = \,\alpha.\]
  • Diagonal entries: \[(A)_{j,j} = 1 + \beta\,(t_0 + V_j), \quad (B)_{j,j} = 1 – \beta\,(t_0 + V_j).\]

Numerical Solution via the Thomas Algorithm

Direct inversion of the \(N\times N\) matrix \(A\) at every time step costs \(O(N^3)\). However, since \(A\) is tri-diagonal, the Thomas algorithm (a specialized Gaussian elimination) solves \(A\,\boldsymbol{\psi} = \mathbf d\) in \(O(N)\) time:

  • Forward sweep: compute modified coefficients \(c’_j\) and right-hand side \(d’_j\).
  • Backward substitution: recover \(\psi_j\) from \(d’_j\).

    We precompute the \(c’_j\) once before time stepping, as the matrix \(A\) does not change.

    Figure: Normalized probability density of the initial Gaussian wave packet (magenta solid) centered at \(x_0=-200\) nm, overlaid with the square potential barrier (black dashed) of width 40 nm and height \(0.85\,\hbar^2k_0^2/(2m)\).

    Python Code

    Below is a complete Python implementation. Comments highlight key steps.

    Wave Packet Simulation 
    
import numpy as np
import matplotlib.pyplot as plt

# Enable interactive plotting
plt.ion()

# Constants
hbar = 1.0545718e-34      # Reduced Planck constant (J·s)
m_e  = 9.10938356e-31     # Electron mass (kg)

# de Broglie wavelength & wave number
wavelength = 5e-9         # meters
k0 = 2 * np.pi / wavelength

# Spatial grid
Lx = 1e-6                 # total length (m)
Nx = 4024
dx = Lx / Nx
x  = np.linspace(-Lx/2, Lx/2, Nx)

# Time grid (for stability)
dt = 0.2e-15              # seconds
Nt = 15000                # number of time steps

# Initial Gaussian wave packet
sigma_x = 60e-9           # width (m)
x0      = -2e-7           # initial center (m)
kx0     = k0
psi0 = np.exp(-(x - x0)**2 / (2*sigma_x**2)) * np.exp(1j * kx0 * x)
psi0 /= np.sqrt(np.sum(np.abs(psi0)**2) * dx)

# Potential: square barrier at center
barrier_width  = 4e-8     # meters
barrier_height = 0.85 * hbar**2 * k0**2 / (2*m_e)
V = np.zeros_like(x)
V[np.abs(x) < barrier_width/2] = barrier_height

# Precompute finite‐difference prefactors
alpha = 1j * hbar * dt / (4 * m_e * dx**2)
beta  = 1j * dt / (2 * hbar)
t0    = hbar**2 / (m_e * dx**2)    # 2*(ħ²/(2m dx²))

# Build Crank–Nicolson tri-diagonal operators A and B
off_diag_A = -alpha
off_diag_B =  alpha              # ← corrected: not alpha.conjugate()
diag_A     = 1 + beta * ( t0 + V )
diag_B     = 1 - beta * ( t0 + V )

# Build arrays for Thomas algorithm
a = np.full(Nx-1, off_diag_A, dtype=complex)  # sub-diagonal of A
b = diag_A.copy()                             # main diagonal of A
c = np.full(Nx-1, off_diag_A, dtype=complex)  # super-diagonal of A

# Precompute modified super-diagonal coefficients (c_prime)
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):
    """Solve A ψ = d for ψ via the Thomas algorithm."""
    # Forward sweep to compute d_prime
    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
    # Back substitution to get ψ
    ψ = np.zeros_like(d, dtype=complex)
    ψ[-1] = d_prime[-1]
    for i in range(Nx-2, -1, -1):
        ψ[i] = d_prime[i] - c_prime[i] * ψ[i+1]
    return ψ

# Initialize wavefunction
psi = psi0.copy()

# Time evolution loop
for t in range(Nt):
    # Build RHS vector d = B @ ψ^n
    d = diag_B * psi
    d[1:]   += off_diag_B * psi[:-1]
    d[:-1]  += off_diag_B * psi[1:]
    
    # Solve for ψ^{n+1}
    psi = thomas_solve(d)
    
    # Dirichlet boundary conditions (ψ=0 at edges)
    psi[0]  = 0+0j
    psi[-1] = 0+0j
    
    # Renormalize to enforce ∑|ψ|² dx = 1
    norm = np.sqrt(np.sum(np.abs(psi)**2) * dx)
    psi /= norm
    
    # Compute transmission probability
    Rho = np.abs(psi)**2
    transmission = np.sum(Rho[x > barrier_width/2]) * dx
    
    # Plot every 200 steps
    if t % 100 == 0:
        plt.clf()
        plt.plot(x*1e9, Rho/np.max(Rho), 'm-', lw=2)
        plt.plot(x*1e9, V/np.max(V), 'k-', lw=2)
        plt.title(f"Step {t}/{Nt}  Trans={transmission:.4f}")
        plt.xlabel("x (nm)")
        plt.ylabel("Probability density")
        plt.pause(0.1)

# Finalize plot
plt.ioff()
plt.show()

    Simulation Parameters

    • Grid: \(L_x = 1\,\mu\text{m}\), \(N_x=4024\) → \(\Delta x \approx 0.25\ \text{nm}\).
    • Time step: \(\Delta t = 0.2\times10^{-15}\,\mathrm{s}\), total steps \(N_t=15{,}000\).
    • Wave packet: initial width \(\sigma_x=60\,\mathrm{nm}\), centered at \(-200\,\mathrm{nm}\) with momentum \(k_0=2\pi/5\,\mathrm{nm}^{-1}\).
    • Barrier: width \(40\,\mathrm{nm}\), height \(0.85\,(\hbar^2k_0^2/2m)\).

    These choices ensure good spatial resolution of the wave packet and stable time stepping under CN’s unconditional stability.

    Plot of the final probability density of a quantum wave packet after tunneling through a potential barrier, showing the transmitted peak beyond the barrier and the reflected component, with horizontal axis labeled position \(x\) and vertical axis labeled \(|\psi(x)|^2\).
    Figure: Final wave packet after tunneling through the potential barrier.

    Results and Discussion

    • Unconditional stability: unlike explicit schemes, CN remains stable for large \(\Delta t\), though accuracy demands \(\Delta t\) small enough to resolve rapid phase oscillations.
    • Reflection and transmission: The plotted probability density clearly shows partial reflection at the barrier and a transmitted tail. The computed transmission coefficient stabilizes to a constant after the wave packet passes the barrier.
    • Performance: Thomas algorithm scales linearly in \(N_x\), making large grids (millions of points) feasible on a modern workstation.

    Conclusion

    The Crank–Nicolson scheme provides a robust, accurate, and stable framework for simulating time-dependent quantum dynamics with arbitrary potentials. By reducing the update at each step to solving a simple tri-diagonal 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.

    Feel free to adapt the code—change the potential, extend to higher dimensions via alternating-direction CN, or couple to time-dependent fields—to explore a wide range of quantum phenomena.

    Leave a Comment

    Your email address will not be published. Required fields are marked *

    Related Posts

    Scroll to Top