Numerical Schrödinger Time Evolution

Numerical Schrödinger propagation evolves a represented state under a represented Hamiltonian. Its central reliability question is whether the computed trajectory approximates the intended initial-value problem at the required times and for the required observables. Norm conservation is necessary for closed Hermitian dynamics, but it is not a measure of phase accuracy.

Helpful background. The time-evolution operator explains unitary dynamics and time ordering; representations explains arrays and inner products.

The semidiscrete evolution equation

In an orthonormal, time-independent basis,

$$i\hbar\dot c(t)=H(t)c(t),\qquad c(t_0)=c_0.$$

A nonorthogonal fixed basis gives $i\hbar S\dot c=Hc$, with norm $c^\dagger Sc$. Orthogonalize once when suitable, or use an integrator that preserves the correct metric.

If $|\psi(t)\rangle=\sum_j c_j(t)|\phi_j(t)\rangle$ uses a moving basis, differentiating both the coefficients and basis functions gives

$$i\hbar S\dot c=(H-i\hbar D)c, \qquad D_{ij}=\langle\phi_i|\dot\phi_j\rangle.$$

The extra matrix accounts for motion of the coordinate system itself. Its identity $\dot S=D+D^\dagger$ ensures conservation of $c^\dagger S c$ for the projected Hermitian dynamics. Omitting $D$ can change the physical state even when the basis only rotates or changes phase.

For a simple coordinate check, let the physical Hamiltonian vanish and use $|\phi_1(t)\rangle=e^{\kappa t}|0\rangle$, $|\phi_2(t)\rangle=|1\rangle$, with real constant $\kappa$ of inverse-time units. The fixed physical state $(|0\rangle+|1\rangle)/\sqrt2$ has

$$S=\operatorname{diag}(e^{2\kappa t},1),\qquad D=\operatorname{diag}(\kappa e^{2\kappa t},0),\qquad c(t)=\frac1{\sqrt2}\begin{pmatrix}e^{-\kappa t}\\1\end{pmatrix}.$$

The moving-basis equation gives $\dot c_1=-\kappa c_1$, keeping $c^\dagger Sc=1$ even though $c^\dagger c=(1+e^{-2\kappa t})/2$ changes. Omitting $D$ incorrectly keeps the coefficients fixed and changes the physical state. This is a coordinate error, not physical amplification or decay.

Returning to a fixed orthonormal basis, Hermitian $H(t)$ gives $d(c^\dagger c)/dt=0$. The energy instead obeys

$$\frac{d}{dt}\langle H\rangle =\left\langle\frac{\partial H}{\partial t}\right\rangle.$$

Thus energy is conserved for a fixed Hamiltonian, but a driven closed system can absorb or release energy. For a time-independent observable $A$,

$$\frac{d\langle A\rangle}{dt} =\frac{i}{\hbar}\langle[H,A]\rangle.$$

These identities provide complementary tests of the propagated state and operator conventions.

Choosing a propagation representation

For a moderate fixed matrix, matrix exponentiation or diagonalization gives an accurate finite-dimensional reference. For large sparse or matrix-free operators, Krylov evolution approximates the exponential action. If kinetic and potential terms are separately easy to exponentiate, split operators exploit that structure.

Crank–Nicolson trades linear solves for stability and norm preservation. Runge–Kutta methods offer flexible error control but generally do not preserve unitarity exactly. For smooth explicitly time-dependent Hamiltonians, Magnus methods incorporate noncommuting time dependence.

The fastest method depends on spectral width, requested time interval, operator-action cost, and accuracy. Compare methods for the same observable tolerance and representation, rather than at the same step size.

Error budgets and trajectory comparisons

Separate spatial or basis error, time-discretization error, exponential-action or linear-solve error, boundary error, and roundoff. Refine one control while keeping the others smaller. A finer spatial grid often increases the largest kinetic eigenvalue, which can force smaller explicit time steps even when low-energy observables have converged.

For normalized pure states, a phase-independent distance is

$$d_{\rm ray}(\psi,\phi) =\min_\alpha\|\psi-e^{i\alpha}\phi\| =\sqrt{2-2|\langle\psi|\phi\rangle|}.$$

Also inspect physical observables and, when interference against a reference matters, the relevant relative phase. Renormalizing each step can conceal norm drift while leaving a wrong trajectory.

A forward-then-backward propagation test checks reversibility of the numerical method, but a reversible approximate method can still be inaccurate. Compare against an analytic solution or a refined trajectory.

A two-state trajectory

For $H=\hbar\Omega\sigma_x/2$ and $c_0=(1,0)^T$,

$$c(t)= \begin{pmatrix} \cos(\Omega t/2)\\ -i\sin(\Omega t/2) \end{pmatrix}, \qquad P_1(t)=\sin^2(\Omega t/2).$$

Norm and energy are constant, yet an incorrect value of $\Omega$ would preserve both. Sampling a complete oscillation and comparing the transition probability detects the phase-frequency error. This benchmark also tests the factor of $\hbar$, angular-frequency units, and multiplication by $i$.

Exercises

If a numerical method preserves norm exactly but accumulates a relative phase error $\delta$ between two equally weighted energy components, what is the fidelity with the exact state?

Solution

After removing the common phase, compare $(1,1)^T/\sqrt2$ with $(1,e^{i\delta})^T/\sqrt2$. The squared overlap is $\cos^2(\delta/2)$, so infidelity is $\delta^2/4+O(\delta^4)$. Exact normalization does not prevent this error.

References

• 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 — algorithm comparison at controlled accuracy.
• Saad, Y., “Analysis of some Krylov subspace approximations to the matrix exponential operator,” SIAM Journal on Numerical Analysis 29, 209–228 (1992), doi:10.1137/0729014 — exponential-action approximation.