Chebyshev Propagation

Chebyshev propagation expands a fixed-Hamiltonian propagator in polynomials whose magnitudes are bounded on the rescaled spectrum. It uses only Hamiltonian–vector products and can achieve high accuracy without storing a Krylov basis. Its essential precondition is a reliable spectral interval: an underestimated bound can make the recurrence unstable.

Helpful background. Matrix exponentials supplies the target action, and Hamiltonian truncation explains why a represented spectrum is bounded.

Spectral scaling and the expansion

Suppose $H=H^\dagger$ has spectrum in $[E_{\min},E_{\max}]$. Set

$$a=\frac{E_{\max}-E_{\min}}2,\qquad b=\frac{E_{\max}+E_{\min}}2,\qquad \widetilde H=\frac{H-bI}{a}.$$

For $a>0$, the spectrum of $\widetilde H$ lies in $[-1,1]$. The Jacobi–Anger expansion gives

$$e^{-iHt/\hbar} =e^{-ibt/\hbar}\left[ J_0(z)I+2\sum_{k=1}^{\infty}(-i)^kJ_k(z)T_k(\widetilde H) \right], \qquad z=\frac{at}{\hbar}.$$

Here $J_k$ are Bessel functions and $T_k$ are Chebyshev polynomials. The scalar phase $e^{-ibt/\hbar}$ restores the energy center. If $a=0$, the Hamiltonian is a scalar on the represented space and propagation is just that phase; do not divide by zero.

Compute vectors

$$v_0=\psi_0,\qquad v_1=\widetilde H\psi_0,\qquad v_{k+1}=2\widetilde H v_k-v_{k-1},$$

and accumulate the weighted sum. Only a few work vectors and the output accumulator are required. Compensated summation or appropriate precision can help when many terms cancel.

Choosing the polynomial degree

Because $\|T_k(\widetilde H)\|_2\leq1$ for Hermitian $\widetilde H$ in the interval, truncating after degree $K$ gives the bound

$$\|U-U_K\|_2 \leq2\sum_{k>K}|J_k(z)|.$$

Choose $K$ from a bound or a carefully evaluated tail, not from one small coefficient: individual Bessel coefficients have zeros before the tail is negligible. The number of significant terms grows roughly with $|z|$, followed by a transition and rapidly decaying tail. Longer times or a wider spectrum therefore require higher degree.

A truncated polynomial is not exactly unitary. Norm defects should decrease with the tail tolerance until roundoff or errors in Hamiltonian actions dominate. A coarse but exactly Hermitian spatial representation does not itself cause norm loss; its physical error requires a separate refinement test. Repeated short intervals and one long polynomial have different roundoff and setup costs; compare at a fixed total accuracy.

Spectral bounds and failure modes

Conservative Gershgorin bounds, analytical kinetic-plus-potential bounds, or validated extremal estimates can supply the interval. Extremal Ritz values from a small subspace are generally inside the true extremal spectrum and are not automatically safe bounds. Add a justified margin or verify the extremal calculation.

Outside $[-1,1]$, $T_k(x)$ grows exponentially with $k$. Even a weakly occupied out-of-range mode can be amplified. Monitor recurrence norms and verify that unexpectedly rapid growth is not caused by a too-narrow interval. The bounded-polynomial argument also relies on Hermiticity; it cannot be transferred unchanged to a non-normal absorbing Hamiltonian.

For a time-dependent Hamiltonian, a frozen-step Chebyshev action must be combined with time-ordering control. Increasing $K$ only solves that frozen action more accurately.

A two-energy benchmark

For $H=\hbar\Omega\sigma_z$, choose $a=\hbar|\Omega|$ and $b=0$. With positive $\Omega$, $\widetilde H=\sigma_z$ and its eigenvalues are $\pm1$. Thus $T_k(1)=1$ and $T_k(-1)=(-1)^k$, and the two accumulated scalar series must reproduce $e^{-i\Omega t}$ and $e^{i\Omega t}$.

An equal superposition tests their relative phase. Repeating with deliberately enlarged bounds should change the required polynomial degree but not the converged state. Repeating with underestimated bounds reveals the instability mechanism.

Exercises

Why can shifting the Hamiltonian reduce work while leaving the relevant spectral half-width unchanged?

Solution

Centering removes the scalar phase from the recurrence and keeps $\widetilde H$ in a symmetric interval. The polynomial degree depends mainly on the half-width $a$, not the absolute energy offset $b$. A poorly chosen interval around zero may be much wider than the centered spectral interval and waste work.

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 — accuracy and efficiency comparisons.
• Tal-Ezer, H., and R. Kosloff, “An accurate and efficient scheme for propagating the time dependent Schrödinger equation,” Journal of Chemical Physics 81, 3967–3971 (1984), doi:10.1063/1.448136 — Chebyshev quantum propagation.