Time-Dependent Schrödinger Equation

Pdes

The time-dependent Schrödinger equation (TDSE)

i ℏ \partial_{t} ψ (r, t) = \hat{H} ψ (r, t), \hat{H} = - \frac{ℏ ^{2}}{2 m} \nabla^{2} + V (r, t)

governs the evolution of a quantum wavefunction under a (possibly time-dependent) Hamiltonian. It is a LINEAR HYPERBOLIC-like PDE (first-order in time, second-order in space, complex-valued), with the distinguishing structural feature that the evolution is UNITARY: it preserves the $L^{2}$ norm of $ψ$ . Any decent numerical scheme should preserve this property too. Two methods dominate: SPLIT-STEP FOURIER (Strang splitting in position/momentum) and CRANK-NICOLSON (implicit, mass-matrix-based). Both are second-order accurate, unconditionally stable, and exactly norm-conserving (to machine precision). This page covers split-step Fourier; Crank-Nicolson is its standard counterpart for problems with non-trivial boundary conditions.

Why TDSE is special

Two structural facts dictate the numerical approach.

(1) Unitarity. The exact solution is $ψ (t) = e^{- i \hat{H} t /ℏ} ψ (0)$ , a unitary evolution operator. The norm $⟨ ψ ∣ ψ ⟩$ is exactly conserved. Numerical schemes that violate this conservation produce wavefunctions whose norm drifts in time — physically meaningless. Methods that EXACTLY preserve unitarity (split-step, Crank-Nicolson) are favored even over higher-order methods that don't.

(2) Operator structure: kinetic + potential. Write $\hat{H} = \hat{T} + \hat{V}$ with $\hat{T} = - ℏ^{2} \nabla^{2} / (2 m)$ the kinetic operator and $\hat{V} = V (r)$ the potential. The kinetic operator is DIAGONAL in momentum space (its Fourier transform is $ℏ^{2} k^{2} / (2 m)$ , a multiplication operator), the potential is DIAGONAL in position space. Each is trivial to exponentiate in its OWN natural basis. They DON'T commute ( $[\hat{T}, \hat{V}] = - i ℏ (\nabla V) \cdot \hat{p} / m \neq = 0$ in general), so splitting incurs an error — but the error vanishes as $Δ t^{3}$ per Strang step.

The split-step Fourier method

Apply Strang splitting to the formal evolution operator:

e^{- i \hat{H} Δ t /ℏ} \approx e^{- i \hat{V} Δ t / (2ℏ)} e^{- i \hat{T} Δ t /ℏ} e^{- i \hat{V} Δ t / (2ℏ)} .

Each step is a sequence of cheap operations on the wavefunction. With $ψ$ sampled on an $N$ -point uniform grid and $ℏ = m = 1$ (atomic units):

Multiply pointwise: $ψ (x) \leftarrow e^{- iV (x) Δ t /2} ψ (x)$ . Half-step in the potential. $O (N)$ .
FFT: $\tilde{ψ} (k) = F [ψ]$ . $O (N lo g N)$ .
Multiply pointwise in momentum: $\tilde{ψ} (k) \leftarrow e^{- i k^{2} Δ t /2} \tilde{ψ} (k)$ . Full step in the kinetic. $O (N)$ .
Inverse FFT: $ψ (x) = F^{- 1} [\tilde{ψ}]$ . $O (N lo g N)$ .
Multiply again pointwise: $ψ (x) \leftarrow e^{- iV (x) Δ t /2} ψ (x)$ . Second half-step. $O (N)$ .

Each step is exactly UNITARY (every operation is multiplication by a complex number of modulus 1 or an FFT, all unitary), so the norm is preserved to machine precision. The accuracy is $O (Δ t^{2})$ in time from Strang and spectral in space from the FFT (the kinetic operator acts EXACTLY on band-limited wavefunctions). For wavefunctions that are smooth on the grid, you get essentially exact spatial discretization.

Two caveats. (1) The method assumes PERIODIC boundary conditions (the FFT requires them). For non-periodic problems, use absorbing boundary layers or the Crank-Nicolson method instead. (2) The momentum grid is determined by the FFT — wave numbers above $π /Δ x$ alias. If the wavefunction develops high-momentum components beyond this limit (e.g. through a sharp tunneling barrier), the simulation aliases unphysically.

Code: harmonic-oscillator test

The harmonic oscillator $V = x^{2} /2$ in atomic units has eigenvalues $E_{n} = n + 1/2$ and a special structure: coherent states (Gaussian wave packets displaced from the origin) oscillate as RIGID BODIES, with the centre executing exact classical motion $⟨ x ⟩ (t) = x_{0} cos t$ . Two tests for the propagator:

EVOLVE THE GROUND STATE: should pick up a global phase $e^{- i E_{0} T}$ over time and return to itself after one period (up to that phase). Check the wavefunction's deviation and its norm.
EVOLVE A COHERENT STATE: should track classical motion. Check $⟨ x ⟩ (t)$ at several phases of the oscillation.

# Split-step Fourier method for the 1D time-dependent Schrödinger equation:
#   i d psi / dt = -1/2 d² psi / dx² + V(x) psi.
#
# Strang splitting in position-momentum representation:
#   psi(t + dt) = exp(-i V dt / 2) F^{-1} exp(-i k² dt / 2) F exp(-i V dt / 2) psi(t)
# where F is the FFT. Each sub-step is diagonal in its natural basis, so
# the algorithm is O(N log N) per time step, unconditionally stable, and
# EXACTLY UNITARY (norm-preserving to machine precision).
#
# Test: SHO ground state phase, and a coherent-state wave packet that
# oscillates classically with period T = 2 pi in V = x²/2 units.

import numpy as np

L  = 20.0                                  # spatial extent
N  = 1024                                  # grid points
x  = np.linspace(-L/2, L/2, N, endpoint=False)
h  = x[1] - x[0]
k  = np.fft.fftfreq(N, d=h) * 2 * np.pi    # momentum-space grid
V  = 0.5 * x**2                            # harmonic-oscillator potential
dt = 0.01

half_V = np.exp(-1j * V * dt / 2)
full_K = np.exp(-1j * 0.5 * k**2 * dt)

def step(psi):
    """One Strang time step."""
    psi = half_V * psi
    psi = np.fft.ifft(full_K * np.fft.fft(psi))
    psi = half_V * psi
    return psi

# ─── Ground-state phase evolution ──────────────────────────────────────
# psi_0(x) = pi^{-1/4} exp(-x²/2),  E_0 = 1/2,  so psi(t) = psi_0 exp(-i t/2).
psi = (1/np.pi**0.25) * np.exp(-x**2 / 2)
T_period = 2 * np.pi
nsteps = int(T_period / dt)
for _ in range(nsteps):
    psi = step(psi)
exact_final = (1/np.pi**0.25) * np.exp(-x**2 / 2) * np.exp(-1j * T_period / 2)
print(f"SHO ground state, one full period T = 2π:")
print(f"  ||psi_final - psi_exact||_inf = {np.max(np.abs(psi - exact_final)):.2e}")
print(f"  <psi|psi> (should be 1)        = {np.trapz(np.abs(psi)**2, x):.10f}")

# ─── Coherent-state wave packet ────────────────────────────────────────
# Initial displacement x0 = 2, zero momentum. Classical oscillation:
# <x>(t) = x0 cos(t). Should track exactly for the SHO (no spreading
# beyond the natural Gaussian width since coherent states are eigenstates
# of the annihilation operator).
print(f"\nCoherent-state wave packet centred at x_0 = 2:")
print(f"  {'t/T':>6s}  {'<x>(t) numerical':>20s}  {'<x>(t) classical':>20s}")
psi0 = (1/np.pi**0.25) * np.exp(-(x - 2.0)**2 / 2)
for frac in [0.0, 0.25, 0.50, 0.75, 1.00]:
    psi_t = psi0.copy()
    n_t = int(frac * T_period / dt)
    for _ in range(n_t):
        psi_t = step(psi_t)
    x_mean = np.trapz(np.abs(psi_t)**2 * x, x)
    classical = 2.0 * np.cos(frac * T_period)
    print(f"  {frac:>6.2f}  {x_mean:>20.6f}  {classical:>20.6f}")

Output:

SHO ground state, one full period T = 2π:
  ||psi_final - psi_exact||_inf = 1.19e-03
  <psi|psi> (should be 1)        = 1.0000000000

Coherent-state wave packet centred at x_0 = 2:
     t/T      <x>(t) numerical     <x>(t) classical
    0.00              2.000000              2.000000
    0.25              0.001580              0.000000
    0.50             -1.999998             -2.000000
    0.75             -0.004739             -0.000000
    1.00              1.999990              2.000000

Three things to read off. (1) Ground-state evolution: error after one full period is $1.2 \times 1 0^{- 3}$ — limited by the $Δ t^{2}$ Strang error accumulating over 628 time steps. Refining $Δ t$ by 10× would drop the error to $\sim 1 0^{- 5}$ (second-order convergence). (2) Norm conservation: $⟨ ψ ∣ ψ ⟩ = 1$ to ten decimal places — exact at machine precision; the global complex phase doesn't leak any norm. (3) Coherent state oscillates classically: $⟨ x ⟩$ at $t / T = 0.25$ and $t / T = 0.75$ is $\sim 1 0^{- 3}$ instead of zero — comparable to the global accuracy. At $t / T = 0.5$ it has flipped to $- 2$ as classical mechanics demands (perfect half-period reflection through the well minimum). At $t / T = 1$ it returns to $2$ .

Crank-Nicolson alternative

For problems with non-periodic boundary conditions (atoms in a box, scattering off a barrier with absorbing boundaries), the standard alternative is Crank-Nicolson:

(I + \frac{i Δ t}{2ℏ} \hat{H}) ψ^{n + 1} = (I - \frac{i Δ t}{2ℏ} \hat{H}) ψ^{n} .

Apply finite differences (or finite elements) to $\hat{H}$ , and solve the tridiagonal linear system at every time step. Crank-Nicolson is also second-order, unconditionally stable, exactly unitary (the operator $(I + i Δ t \hat{H} /2)^{- 1} (I - i Δ t \hat{H} /2)$ is the Cayley transform of $- i Δ t \hat{H} /2$ , manifestly unitary), and handles arbitrary boundary conditions naturally. Cost per step is $O (N)$ for a tridiagonal system in 1D, $O (N^{3/2})$ for sparse 2D systems via banded solvers — slightly more expensive than split-step Fourier per step, but no periodic-domain constraint.

Applications

Quantum dynamics in molecular systems: model nuclear motion on a Born-Oppenheimer potential surface; tunneling, photo-dissociation, vibrational dynamics. The standard tool of quantum-chemical wave-packet dynamics.
Strong-field physics: above-threshold ionisation, high-harmonic generation. The TDSE with a time-dependent laser-pulse vector potential needs split-step or Crank-Nicolson on a fine grid.
Wave-packet dynamics in 2D semiconductors: time-resolved electron dynamics in graphene, surface states, topological materials. Split-step Fourier in periodic 2D is a workhorse.
Schrödinger-Poisson: the TDSE coupled to its own Hartree potential, a mean-field approximation to many-body quantum dynamics. Time-dependent DFT is the next step up (cf. Liouville-Lanczos for the linear-response version).
Real-time TDDFT: full ab-initio quantum dynamics of electrons. Same numerical framework as the simple TDSE, just with much more complicated $V (r, t)$ coming from the Kohn-Sham equations.

Connections to the rest of the site

The time-dependent Schrödinger equation sits at the intersection of several themes already on the site. The quantum physics section covers the conceptual framework. The harmonic oscillator (finite difference) page covers the time-INDEPENDENT eigenvalue problem on the same potential we used above. The WKB page covers the semiclassical limit of the same equation. The Liouville-Lanczos page covers the linear-response framework for excitation spectra. Each tackles a different facet of the same underlying equation; together they form a complete picture of how the Schrödinger equation is solved in practice.

Operator splitting — the technique underlying split-step Fourier.
Spectral methods — the FFT spatial discretization gives spectral accuracy in space.
The three classical equations — TDSE behaves like the wave equation (first-order in time but complex, so reflects as a phase rotation).
WKB — the $ℏ \to 0$ semiclassical limit of the equation we're solving here.
Harmonic oscillator (FD) — the time-independent counterpart on the same potential.