A beginner friendly introduction to Chebyshev Propagator

1. When is a Chebyshev propagator needed?

1.1 Time evolution as an operator exponential

In quantum mechanics (setting $\hbar = 1$), a time-independent Hamiltonian $\hat{H}$ generates dynamics via the Schrödinger equation

\[ i\frac{d}{dt}\,|\psi(t)\rangle = \hat{H}\,|\psi(t)\rangle. \tag{1} \]

The formal solution is

\[ |\psi(t)\rangle = e^{-i\hat{H}t}|\psi(0)\rangle. \tag{2} \]

Thus, the central computational task is:

Task: Given $\hat{H}$ and $|\psi(0)\rangle$, compute $|\psi(t)\rangle = e^{-i\hat{H}t}|\psi(0)\rangle$.

1.2 Why we do not want to diagonalize

If the Hilbert space dimension is large, diagonalizing $\hat{H}$ or forming the dense matrix $e^{-i\hat{H}t}$ is impractical. However, we only need the action of the propagator on a vector:

\[ |\psi(t)\rangle = \big(e^{-i\hat{H}t}\big)\,|\psi(0)\rangle. \tag{3} \]

Chebyshev propagation achieves this by approximating the operator exponential by a polynomial series in $\hat{H}$ (after rescaling), and evaluating the series using a stable three-term recurrence that requires only repeated matrix–vector multiplications.

1.3 Layout

We proceed in the following connected steps:

1. Define Chebyshev polynomials $T_n(x)$, prove boundedness on $[-1,1]$, and derive their recurrence.
2. Explain why the Hamiltonian must be rescaled so its spectrum lies in $[-1,1]$.
3. Derive the scalar Chebyshev expansion of $e^{-i\tau x}$ for $x\in[-1,1]$: \[ e^{-i\tau x} = J_0(\tau) + 2\sum_{n\ge 1}(-i)^n J_n(\tau) T_n(x). \]
4. Lift the scalar identity to an operator identity $x\mapsto \hat{H}'$.
5. Derive the Chebyshev propagation algorithm (vector recurrence).
6. Derive a truncation error bound and explain why $J_n(\tau)$ becomes small for $n\gtrsim \tau$.
7. Work through explicit $2\times 2$ examples to demonstrate how the propagator is used.

3. Why we rescale the Hamiltonian to $[-1,1]$

3.1: Eigenvalues determine operator functions

Let $\hat{H}$ be Hermitian. It has an orthonormal eigenbasis $\{|E_k\rangle\}$ with real eigenvalues $\{E_k\}$:

\[ \hat{H}|E_k\rangle = E_k|E_k\rangle,\qquad E_k\in\mathbb{R}. \tag{13} \]

The propagator acts on eigenstates as

\[ e^{-i\hat{H}t}|E_k\rangle = e^{-iE_k t}|E_k\rangle. \tag{14} \]

So approximating $e^{-i\hat{H}t}$ reduces to approximating $e^{-iEt}$ over the spectrum.

3.2: Why Chebyshev polynomials demand $x\in[-1,1]$

Chebyshev polynomials are bounded on $[-1,1]$ by (10). Outside that interval, $T_n(x)$ grows rapidly with $n$ (in fact exponentially). Thus, for stability of $T_n(\hat{H}')$ at large $n$, we want $\mathrm{spec}(\hat{H}')\subset[-1,1]$.

3.3: Rescaling the spectrum

Assume bounds $E_{\min}\le E_k\le E_{\max}$ for all eigenvalues:

\[ E_{\min}\le E_k \le E_{\max}. \tag{15} \]

Define

\[ a=\frac{E_{\max}+E_{\min}}{2},\qquad b=\frac{E_{\max}-E_{\min}}{2}>0, \tag{16} \]

and the shifted/rescaled Hamiltonian

\[ \hat{H}'=\frac{\hat{H}-a\hat{I}}{b}. \tag{17} \]

If $\hat{H}|E_k\rangle=E_k|E_k\rangle$, then

\[ \hat{H}'|E_k\rangle=\frac{E_k-a}{b}|E_k\rangle \equiv E_k'|E_k\rangle. \tag{18} \]

Check endpoints: $E_{\max}\mapsto +1$ and $E_{\min}\mapsto -1$, hence

\[ E_k'\in[-1,1]. \tag{19} \]

3.4: Factorizing the propagator

Rewrite $\hat{H}=a\hat{I}+b\hat{H}'$. Then

\[ e^{-i\hat{H}t}=e^{-i(a\hat{I}+b\hat{H}')t}. \tag{20} \]

Since $\hat{I}$ commutes with $\hat{H}'$, we can factor:

\[ e^{-i\hat{H}t}=e^{-iat}\,e^{-i(bt)\hat{H}'}. \tag{21} \]

Define the scaled time

\[ \tau\equiv bt. \tag{22} \]

Now the numerical problem is computing $e^{-i\tau \hat{H}'}$ with spectrum inside $[-1,1]$.

5. Lifting the scalar identity to an operator identity

5.1: Eigen-decomposition of $\hat{H}'$

Since $\hat{H}'$ is Hermitian, it has eigenpairs $\{E_k',|E_k\rangle\}$ with $E_k'\in[-1,1]$:

\[ \hat{H}'|E_k\rangle = E_k'|E_k\rangle. \tag{40} \]

Completeness:

\[ \hat{I} = \sum_k |E_k\rangle\langle E_k|. \tag{41} \]

5.2: Action of $e^{-i\tau\hat{H}'}$ on eigenstates

\[ e^{-i\tau\hat{H}'}|E_k\rangle = e^{-i\tau E_k'}|E_k\rangle. \tag{42} \]

5.3: Apply the scalar identity to each eigenvalue

Since $E_k'\in[-1,1]$, apply (39) with $x=E_k'$:

\[ e^{-i\tau E_k'} = J_0(\tau)+2\sum_{n=1}^{\infty}(-i)^n J_n(\tau) T_n(E_k'). \tag{43} \]

Insert into (42):

\[ e^{-i\tau\hat{H}'}|E_k\rangle = \left[ J_0(\tau)+2\sum_{n=1}^{\infty}(-i)^n J_n(\tau) T_n(E_k') \right]|E_k\rangle. \tag{44} \]

5.4: Define $T_n(\hat{H}')$ via spectral calculus

Define the operator polynomial

\[ T_n(\hat{H}') \equiv \sum_k T_n(E_k')\,|E_k\rangle\langle E_k|. \tag{45} \]

Then

\[ T_n(\hat{H}')|E_k\rangle = T_n(E_k')|E_k\rangle. \tag{46} \]

So (44) can be rewritten as

\[ e^{-i\tau\hat{H}'}|E_k\rangle = \left[ J_0(\tau)\hat{I}+2\sum_{n=1}^{\infty}(-i)^n J_n(\tau)\,T_n(\hat{H}') \right]|E_k\rangle. \tag{47} \]

Since the eigenstates form a basis, the operator identity follows:

\[ e^{-i\tau\hat{H}'} = J_0(\tau)\hat{I}+2\sum_{n=1}^{\infty}(-i)^n J_n(\tau)\,T_n(\hat{H}'). \tag{48} \]

Finally, restore the original Hamiltonian using (21) and $\tau=bt$:

\[ e^{-i\hat{H}t} = e^{-iat}\left[ J_0(bt)\hat{I}+2\sum_{n=1}^{\infty}(-i)^n J_n(bt)\,T_n(\hat{H}') \right]. \tag{49} \]

7. Truncation error and why $J_n(\tau)$ gets small for $n\gtrsim\tau$

7.1: Define the remainder operator

Let

\[ \hat{U}'(\tau)\equiv e^{-i\tau\hat{H}'}. \tag{62} \]

Define the truncated approximation

\[ \hat{U}'_N(\tau) \equiv J_0(\tau)\hat{I} + 2\sum_{n=1}^{N}(-i)^n J_n(\tau)\,T_n(\hat{H}'). \tag{63} \]

The remainder is

\[ \hat{R}_N(\tau)\equiv \hat{U}'(\tau)-\hat{U}'_N(\tau) = 2\sum_{n=N+1}^{\infty}(-i)^n J_n(\tau)\,T_n(\hat{H}'). \tag{64} \]

7.2: Prove $\|T_n(\hat{H}')\|\le 1$

In the eigenbasis of $\hat{H}'$, the eigenvalues of $T_n(\hat{H}')$ are $T_n(E_k')$. Since $E_k'\in[-1,1]$, boundedness gives $|T_n(E_k')|\le 1$. Therefore

\[ \|T_n(\hat{H}')\|\le 1. \tag{65} \]

7.3: Derive the error bound

Take norms in (64) and apply triangle inequality:

\[ \|\hat{R}_N(\tau)\| \le 2\sum_{n=N+1}^{\infty}|J_n(\tau)|\,\|T_n(\hat{H}')\| \tag{66} \] \[ \le 2\sum_{n=N+1}^{\infty}|J_n(\tau)|. \tag{67} \]

Thus

\[ \|\hat{U}'(\tau)-\hat{U}'_N(\tau)\| \le 2\sum_{n=N+1}^{\infty}|J_n(\tau)|. \tag{68} \]

Applying to a normalized state $\||\psi(0)\rangle\|=1$:

\[ \||\psi_{\mathrm{exact}}(t)\rangle-|\psi_N(t)\rangle\| \le 2\sum_{n=N+1}^{\infty}|J_n(\tau)|. \tag{69} \]

7.4: Why the Bessel tail becomes small (intuition)

Recall $(-i)^nJ_n(\tau)$ appears as a Fourier coefficient of $e^{-i\tau\cos\theta}$:

\[ (-i)^n J_n(\tau) = \frac{1}{\pi}\int_{0}^{\pi} e^{-i\tau\cos\theta}\cos(n\theta)\,d\theta. \tag{70} \]

For large $n$, $\cos(n\theta)$ oscillates rapidly, producing cancellations. Consider the oscillatory phase

\[ \Phi(\theta)=\tau\cos\theta - n\theta. \tag{71} \]

When $n\gg\tau$, the $-n\theta$ term dominates the $\theta$-variation, so $\Phi(\theta)$ changes rapidly and cancellations are strong. Hence $J_n(\tau)$ becomes small for $n\gtrsim\tau$.

A practical rule is

\[ N \approx \tau + (\text{a modest safety margin}),\qquad \tau=bt. \tag{72} \]

9. Example 2: a nontrivial rescaling example

Let

\[ \hat{H}= \begin{pmatrix} 2 & 1\\ 1 & 4 \end{pmatrix}, \qquad |\psi(0)\rangle= \begin{pmatrix} 1\\ 0 \end{pmatrix}. \tag{82} \]

9.1: Compute eigenvalues exactly

Characteristic polynomial:

\[ \det(\hat{H}-\lambda\hat{I})= \begin{vmatrix} 2-\lambda & 1\\ 1 & 4-\lambda \end{vmatrix} = (2-\lambda)(4-\lambda)-1. \tag{83} \] \[ = \lambda^2 - 6\lambda + 7. \tag{84} \]

Solve $\lambda^2-6\lambda+7=0$:

\[ \lambda_{\pm}=\frac{6\pm\sqrt{36-28}}{2} = 3\pm\sqrt{2}. \tag{85} \]

Thus $E_{\max}=3+\sqrt{2}$, $E_{\min}=3-\sqrt{2}$.

9.2: Compute $a$, $b$ and $\hat{H}'$

\[ a=\frac{E_{\max}+E_{\min}}{2}=3,\qquad b=\frac{E_{\max}-E_{\min}}{2}=\sqrt{2}. \tag{86} \] \[ \hat{H}'=\frac{\hat{H}-a\hat{I}}{b} = \frac{1}{\sqrt{2}} \left[ \begin{pmatrix} 2 & 1\\ 1 & 4 \end{pmatrix} - 3\begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \right] = \frac{1}{\sqrt{2}} \begin{pmatrix} -1 & 1\\ 1 & 1 \end{pmatrix}. \tag{87} \]

Scaled time is $\tau=bt=\sqrt{2}\,t$, and

\[ e^{-i\hat{H}t}=e^{-i3t}\,e^{-i\tau\hat{H}'}. \tag{88} \]

9.3: Compute Chebyshev vectors explicitly

\[ |\phi_0\rangle = \begin{pmatrix} 1\\ 0 \end{pmatrix}. \tag{89} \] \[ |\phi_1\rangle=\hat{H}'|\phi_0\rangle= \frac{1}{\sqrt{2}} \begin{pmatrix} -1 & 1\\ 1 & 1 \end{pmatrix} \begin{pmatrix} 1\\ 0 \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} -1\\ 1 \end{pmatrix}. \tag{90} \]

Now compute $\hat{H}'|\phi_1\rangle$ carefully:

\[ \hat{H}'|\phi_1\rangle = \frac{1}{\sqrt{2}} \begin{pmatrix} -1 & 1\\ 1 & 1 \end{pmatrix} \cdot \frac{1}{\sqrt{2}} \begin{pmatrix} -1\\ 1 \end{pmatrix} = \frac{1}{2} \begin{pmatrix} (-1)(-1) + 1\cdot 1\\ 1\cdot(-1) + 1\cdot 1 \end{pmatrix} = \frac{1}{2} \begin{pmatrix} 2\\ 0 \end{pmatrix} = \begin{pmatrix} 1\\ 0 \end{pmatrix} = |\phi_0\rangle. \tag{92} \]

Then the recurrence gives

\[ |\phi_2\rangle = 2\hat{H}'|\phi_1\rangle - |\phi_0\rangle = 2|\phi_0\rangle-|\phi_0\rangle = |\phi_0\rangle. \tag{93} \]

Similarly one obtains the same even/odd pattern:

\[ |\phi_{2m}\rangle = |\phi_0\rangle,\qquad |\phi_{2m+1}\rangle = |\phi_1\rangle. \tag{94} \]

9.4: Assemble the propagated state

\[ |\psi(t)\rangle \approx e^{-i3t}\left[ J_0(\tau)|\phi_0\rangle + 2\sum_{n=1}^{N}(-i)^n J_n(\tau)|\phi_n\rangle \right],\qquad \tau=\sqrt{2}\,t. \tag{95} \]

Separating even/odd terms yields a two-vector form:

\[ |\psi(t)\rangle \approx e^{-i3t}\Big[A_N(\tau)|\phi_0\rangle + B_N(\tau)|\phi_1\rangle\Big], \tag{96} \] \[ A_N(\tau)\equiv J_0(\tau) + 2\sum_{m=1}^{\lfloor N/2\rfloor}(-1)^m J_{2m}(\tau), \tag{97} \] \[ B_N(\tau)\equiv -2i\sum_{m=0}^{\lfloor (N-1)/2\rfloor}(-1)^m J_{2m+1}(\tau). \tag{98} \]

In the limit $N\to\infty$, $A_N(\tau)\to\cos\tau$ and $B_N(\tau)\to -i\sin\tau$, so

\[ |\psi(t)\rangle = e^{-i3t}\left[\cos(\sqrt{2}\,t)\,|\phi_0\rangle - i\sin(\sqrt{2}\,t)\,|\phi_1\rangle\right]. \tag{99} \]