Time-Stepping Methods

Time-stepping methods approximate the time-dependent Schrödinger equation by advancing a state through small time increments. After discretization in space or in a finite basis, the problem becomes a large system of ordinary differential equations:

$$\frac{d\psi}{dt} = - \frac{i}{\hbar} H(t)\psi(t).$$

The exact closed-system evolution is norm-preserving. A numerical time step that ignores this structure can look accurate for a short time and still accumulate unphysical norm growth, damping, phase error, or energy drift.

This page compares basic explicit and implicit schemes, explains stability, and highlights unitarity diagnostics. It does not own the general theory of time evolution; see Schrödinger Equation and Unitary Time Evolution for the physics.

Exact One-Step Benchmark

If $H$ is time independent, the exact solution over one step is

$$\psi(t+\Delta t) = \exp\left( - \frac{i}{\hbar} H\Delta t \right) \psi(t).$$

For Hermitian $H$, this operator is unitary:

$$U^\dagger U = I.$$

Thus the norm satisfies

$$\lVert \psi(t+\Delta t)\rVert = \lVert \psi(t)\rVert.$$

Every time-stepping method for closed-system quantum mechanics should be judged against this benchmark. A method may be stable, accurate, cheap, or exactly unitary, but not every method has all of those properties.

Step-Size Notation

Let

$$t_n=t_0+n\Delta t, \qquad \psi_n\approx\psi(t_n).$$

A one-step method advances

$$\psi_n \mapsto \psi_{n+1}.$$

For a time-independent finite Hamiltonian, many methods can be analyzed mode by mode. If $H\phi=E\phi$, then the exact amplification factor over one step is

$$\exp(-iz), \qquad z=\frac{E\Delta t}{\hbar}.$$

Because $\lvert\exp(-iz)\rvert=1$ for real $z$, any numerical amplification factor with systematic modulus different from $1$ introduces artificial damping or growth.

Explicit Euler

The simplest explicit method applies the right-hand side at the current time:

$$\psi_{n+1} = \psi_n - \frac{i\Delta t}{\hbar} H\psi_n.$$

For an energy eigenmode, the numerical amplification factor is

$$g_E = 1-iz.$$

Its squared modulus is

$$\lvert g_E\rvert^2 = 1+z^2.$$

Therefore explicit Euler grows the norm for every nonzero real energy mode. It is not a good real-time Schrödinger integrator except as a warning example or a very short-time local truncation exercise.

Implicit Euler

Implicit Euler evaluates the right-hand side at the next time:

$$\psi_{n+1} = \psi_n - \frac{i\Delta t}{\hbar} H\psi_{n+1}.$$

Equivalently,

$$\left( I+ \frac{i\Delta t}{\hbar} H \right) \psi_{n+1} = \psi_n.$$

For an energy eigenmode,

$$g_I = \frac{1}{1+iz}, \qquad \lvert g_I\rvert^2 = \frac{1}{1+z^2}.$$

This damps the norm. It can be stable for dissipative differential equations, but damping is unphysical for closed real-time quantum dynamics. The method is useful conceptually because it shows that stability alone is not the same as unitary accuracy.

Crank–Nicolson

The Crank–Nicolson method averages the right-hand side between the old and new times:

$$\left( I+ \frac{i\Delta t}{2\hbar} H \right) \psi_{n+1} = \left( I- \frac{i\Delta t}{2\hbar} H \right) \psi_n.$$

For an energy eigenmode,

$$g_{CN} = \frac{1-iz/2}{1+iz/2}.$$

For real $z$,

$$\lvert g_{CN}\rvert = 1.$$

Thus, for time-independent Hermitian $H$ and exact linear solves, Crank–Nicolson is exactly norm-preserving. It is also second-order accurate in $\Delta t$.

The price is an implicit solve at every step. If the linear solve is loose, the numerical evolution will not be exactly unitary even if the formula is unitary in exact arithmetic.

Runge–Kutta Methods

Explicit Runge–Kutta methods evaluate the right-hand side several times inside a step. The familiar fourth-order method is often accurate for small finite systems over moderate times, but it is not exactly unitary.

For a linear autonomous equation $y'=Ay$, a Runge–Kutta method has a stability function $R$:

$$y_{n+1} = R(\Delta t A)y_n.$$

For fourth-order Runge–Kutta,

$$R(z) = 1+z+\frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}.$$

For quantum dynamics, $z$ lies near the imaginary axis when $A=-iH/\hbar$. Good behavior requires $\lvert R(z)\rvert$ to stay close to $1$ over the relevant spectral range. High-energy grid modes can violate this condition even when low-energy observables appear smooth.

Adaptive Runge–Kutta methods are useful for general ODEs and time-dependent finite systems, but their error controllers usually track local truncation error rather than exact quantum unitarity. The general numerical ODE background is ODE Solvers.

Stability and Spectral Range

After spatial discretization, the largest represented energy scale controls the time step. A finite-difference kinetic-energy matrix has high-energy grid modes near the cutoff. Even if the physical state is low energy, numerical noise or rough potentials can excite those modes.

For a method with stability function $R$, inspect

$$\lvert R(-iE\Delta t/\hbar)\rvert$$

over the energies $E$ represented by the finite Hamiltonian. If the method grows high-energy modes, a simulation can fail long before the low-energy physics has a chance to converge.

This is why time-step convergence should be checked together with grid or basis convergence. See Discretization, PDE Solvers, Convergence Tests, and Conditioning and Stability.

Time-Dependent Hamiltonians

When $H(t)$ depends on time, exact evolution is a time-ordered exponential. The order of Hamiltonians at different times matters when

$$[H(t_1),H(t_2)] \ne 0.$$

A simple midpoint exponential approximation is

$$\psi_{n+1} \approx \exp\left( - \frac{i}{\hbar} H\left(t_n+\frac{\Delta t}{2}\right)\Delta t \right) \psi_n.$$

This captures some second-order behavior, but noncommutativity creates additional errors. The conceptual background is Time Ordering.

For rapidly driven systems, compare time-step refinement with physical timescales such as drive periods, pulse widths, avoided-crossing times, and the inverse spectral bandwidth.

Unitarity Diagnostics

For a closed system with Hermitian $H$, monitor

$$N_n = \langle \psi_n,\psi_n\rangle.$$

The deviation

$$N_n-N_0$$

is a basic diagnostic. For time-independent $H$, also monitor energy drift:

$$E_n = \langle \psi_n,H\psi_n\rangle.$$

Norm conservation alone is not sufficient. A method can preserve norm while accumulating phase error or dispersing a wave packet incorrectly. Energy, symmetry quantum numbers, expectation values, and benchmark solutions all give additional checks.

Renormalizing $\psi_n$ after every step can hide a bad integrator. It may be useful as an emergency diagnostic, but it should not be mistaken for a unitary method.

Sparse and Matrix-Free Time Steps

Large quantum dynamics calculations rarely form dense evolution matrices. A step may use:

• sparse matrix-vector products;
• sparse linear solves;
• matrix-free stencil actions;
• basis transformations;
• Krylov approximations to exponential actions.

The storage background is Sparse Matrices. If an implicit method requires solving a linear system, the solve tolerance becomes part of the time-stepping error budget.

Choosing a Method

Situation	Common first choice	Main caution
small dense finite system	exact exponential or high-order ODE solver	phase and long-time error
grid Hamiltonian, closed system	Crank–Nicolson or exponential-action method	solve tolerance and boundary effects
smooth wave packet with Fourier grid	transform-based method	aliasing and grid conventions
time-dependent Hamiltonian	midpoint, Magnus-type, or adaptive ODE method	time ordering and noncommutation
dissipative or open system	problem-specific master-equation integrator	trace, positivity, and stiffness

This table is only a starting point. The correct method depends on the observable, time scale, spectral bandwidth, desired accuracy, and conservation laws.

Common Mistakes

• Using explicit Euler for real-time Schrödinger evolution.
• Calling a method stable without checking norm and phase behavior on imaginary-axis modes.
• Choosing $\Delta t$ from the physical frequency of interest while ignoring high-energy grid modes.
• Treating solver tolerance in an implicit method as unrelated to time-step accuracy.
• Renormalizing after every step and assuming the dynamics became unitary.
• Comparing wavefunctions at long times without separating phase error from shape error.
• Checking time-step convergence while holding an unconverged spatial grid fixed.
• Ignoring time ordering for noncommuting time-dependent Hamiltonians.

Cross-Links

• Schrödinger Equation
• Unitary Time Evolution
• Time-Evolution Operator
• Time Ordering
• Matrix Functions and Exponentials
• Matrix Exponentials Numerically
• Unitary Operators
• ODE Solvers
• PDE Solvers
• Convergence Tests
• Ordinary Differential Equations
• Discretization
• Finite Difference Methods
• Sparse Matrices
• Fast Fourier Transform
• Conditioning and Stability
• Floating-Point Arithmetic
• Time-Dependent Hamiltonian Notebook

References

• E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration, 2nd ed., Springer, 2006.
• E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993.
• R. Kosloff, “Time-dependent quantum-mechanical methods for molecular dynamics”, Journal of Physical Chemistry 92, 2087-2100, 1988.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• T. Pang, An Introduction to Computational Physics, 2nd ed., Cambridge University Press, 2006.
• D. J. Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective, University Science Books, 2007.

Exercises

1. For an energy eigenmode with $z=E\Delta t/\hbar$, show that explicit Euler grows the norm.

Solution

The explicit Euler amplification factor is

$$g_E=1-iz.$$

For real $z$,

$$\lvert g_E\rvert^2 = (1-iz)(1+iz) = 1+z^2.$$

This is greater than $1$ for any nonzero $z$, so repeated steps produce artificial norm growth.

2. Show that the Crank–Nicolson scalar amplification factor has unit modulus for real $z$.

Solution

The factor is

$$g_{CN} = \frac{1-iz/2}{1+iz/2}.$$

The numerator and denominator are complex conjugates up to order, so

$$\lvert g_{CN}\rvert^2 = \frac{1+z^2/4}{1+z^2/4} = 1.$$

3. Why can a time step that resolves the physical oscillation frequency still fail on a fine spatial grid?

Solution

The finite Hamiltonian contains high-energy grid modes near the spatial cutoff. Stability and phase accuracy depend on the whole represented spectral range, not only on the physical frequency one intends to observe. If the time step gives poor amplification for high-energy modes, numerical noise or coupling through a rough potential can make those modes contaminate the calculation.

4. Why is renormalizing the wavefunction after every step not the same as using a unitary method?

Solution

Renormalization fixes only the total norm. It does not correct relative phases, dispersion errors, wrong mode amplitudes, energy drift, broken symmetries, or nonunitary distortion of the state before rescaling. A genuinely unitary step preserves all inner products, not just the norm of one vector after manual adjustment.