Runge–Kutta Methods

Runge–Kutta (RK) methods solve the Schrödinger equation as a system of ordinary differential equations. They are convenient for time-dependent drives and coupled auxiliary equations, but generic explicit RK updates are not unitary. Their local error estimator and stability region must be interpreted in the energy scales of the represented Hamiltonian.

Helpful background. The numerical Schrödinger equation supplies norm and energy identities.

Classical fourth-order stepping

Write $\dot y=f(t,y)=-iH(t)y/\hbar$. Classical RK4 evaluates

$$\begin{aligned} k_1&=f(t,y),\\ k_2&=f(t+h/2,y+hk_1/2),\\ k_3&=f(t+h/2,y+hk_2/2),\\ k_4&=f(t+h,y+hk_3), \end{aligned}$$

and updates $y_{\rm new}=y+h(k_1+2k_2+2k_3+k_4)/6$. For sufficiently smooth problems, the local error is $O(h^5)$ and the global error at fixed final time is $O(h^4)$.

Every stage must use the Hamiltonian at its own stage time. Reusing $H(t)$ throughout a driven step changes the method and generally loses the advertised order. For a complex state, use a solver that supports complex arithmetic or split real and imaginary components consistently.

Stability on quantum eigenmodes

For a fixed eigenmode with $\dot y=\lambda y$, RK4 has stability polynomial

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

Quantum eigenvalues give $\lambda=-iE/\hbar$. For real $x$,

$$|R(ix)|^2=1-\frac{x^6}{72}+\frac{x^8}{576}.$$

Consequently the imaginary-axis stable segment is $|x|\leq2\sqrt2$. Within it, the norm generally drifts downward except at special points; outside it, repeated steps can grow. This is a stability restriction, not an accuracy target. Good phase accuracy usually requires a much smaller step.

A spatial grid can have high kinetic energies proportional to $\Delta x^{-2}$, so refining space can tighten the explicit RK step restriction. A weakly populated high-energy mode may still be excited by numerical noise. Removing an irrelevant scalar energy offset can reduce the range seen by the solver, with its global phase accounted for separately.

Adaptive local tolerances

Embedded RK pairs estimate error by comparing two formulas of different order. A common component scale is $\mathrm{atol}+\mathrm{rtol}\max(|y_n|,|y_{n+1}|)$. In a high-dimensional quantum state, many physically relevant components may be small, so absolute tolerances and norm conventions matter. A per-component threshold is not directly a bound on an observable or total state error.

Specify output sampling separately from internal steps. A solver may take long internal steps between densely requested output times, using interpolation. Force step boundaries at discontinuous pulse changes; a smooth-problem estimator can otherwise miss the discontinuity.

Renormalization changes the numerical map. If used, retain the pre-normalization norm defect and verify phase and observables independently. Renormalizing every step does not turn RK4 into a fourth-order unitary propagator with the same error constants.

A controlled convergence experiment

For $H=\hbar\Omega\sigma_x/2$, propagate $|0\rangle$ for a fixed complete fraction of a Rabi cycle and compare with the analytic state. Use $h,h/2,h/4$ with all other settings fixed. Once the asymptotic regime is reached, global state errors should decrease by about $16$ for RK4. Measure norm drift and phase-sensitive error separately.

A time-dependent noncommuting Hamiltonian is a stronger second test because it exercises stage-time evaluation. Compare it with a refined Magnus or dense short-step exponential calculation.

Exercises

At $x=2$, evaluate the squared RK4 amplification magnitude and explain the consequence over repeated steps.

Solution

$|R(2i)|^2=1-64/72+256/576=5/9$. Each step multiplies the norm by $\sqrt5/3$ and the squared norm by $5/9$. The mode is bounded but strongly damped, even though the exact quantum evolution has neither damping nor growth. Operating inside the stability region is far weaker than achieving accurate evolution.

References

• Hairer, E., S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993 — RK order, stability, and adaptive error control.
• Leforestier, C., et al., “A comparison of different propagation schemes for the time dependent Schrödinger equation,” Journal of Computational Physics 94, 59–80 (1991), doi:10.1016/0021-9991(91)90137-A — quantum propagation accuracy and comparison criteria.