advanced30 min studyLesson 934 in Physics

Crank–Nicolson unitarity and solver tolerance

Analyze a unitary rational propagator and the errors introduced by phase distortion and incomplete linear solves.

On this page

What you’ll understand

  • Derive the Cayley update for a finite Hermitian Hamiltonian.
  • Compute its eigenmode phase and distinguish unitarity from temporal accuracy.
  • Relate a linear-solve residual to a one-step state error.

An implicit method can preserve the norm of a closed finite quantum system while avoiding the growth of explicit Euler propagation. Crank–Nicolson achieves this through a rational approximation to the exponential, often called a Cayley transform. Its useful structure has precise conditions: the Hamiltonian must be Hermitian in the stated inner product, and the implicit system must actually be solved. Even then, finite-timestep phase error remains and can accumulate into observable mistakes.

Derive the rational update from midpoint averaging

For a time-independent finite Hamiltonian, average the right-hand side of the Schrödinger equation between old and new states. Rearranging gives one matrix multiplying the unknown new state and its adjoint multiplying the old state. Since both are polynomials in the same Hermitian H, they commute. Their ratio is unitary in exact arithmetic. Implementations should solve this linear system directly rather than forming a dense matrix inverse merely because inverse notation is convenient.

The norm statement uses the ordinary inner product for an orthonormal finite basis. With an overlap matrix, the update and adjoint relationships must be formulated in that metric. An absorbing non-Hermitian Hamiltonian or a dissipative density-matrix generator belongs to a different analysis. Calling every midpoint-looking update unitary ignores the operator and state-space assumptions that make this particular proof work.

Acn+1=Bcn,A=I+iΔt2ℏH,B=I−iΔt2ℏH=A†,UCN=A−1BA c_{n+1}=B c_n,\qquad A=I+\frac{i\Delta t}{2\hbar}H,\quad B=I-\frac{i\Delta t}{2\hbar}H=A^\dagger,\qquad U_{\rm CN}=A^{-1}B
The matrices are dimensionless. For the declared finite Hermitian H and exact solve, U-CN is unitary and is a second-order approximation over fixed time.↔ Scroll wide equations horizontally.

Measure the eigenphase rather than only the norm

On an eigenvector with real energy E, the update multiplies by a complex number of unit magnitude. Its phase is minus twice the arctangent of E times timestep over twice hbar. The exact phase is instead minus E times timestep over hbar. They agree at small dimensionless timestep through the expected order, but a large step can badly distort relative phases between energy levels even while the norm remains exactly one.

For constant H, the rational propagator also commutes with H and preserves its finite expectation value. Thus norm and energy can both pass while a superposition rotates at the wrong relative rate. For a driven Hamiltonian, using a properly time-centered Hermitian matrix retains per-step unitarity and can give second-order temporal accuracy under appropriate smoothness. Replacing the true time ordering by a poorly sampled midpoint still requires independent step-size checks.

uE=1−iEΔt/(2ℏ)1+iEΔt/(2ℏ)=e−2iarctan⁡[EΔt/(2ℏ)],∣uE∣=1u_E=\frac{1-iE\Delta t/(2\hbar)}{1+iE\Delta t/(2\hbar)}=e^{-2i\arctan[E\Delta t/(2\hbar)]},\qquad |u_E|=1
The eigenmode factor has unit modulus, but its arctangent phase is not generally the exact exponential phase. A global energy offset can also affect finite-step relative phase errors.↔ Scroll wide equations horizontally.

Include iterative-solver error in the propagation budget

Let the linear residual be b minus A times the computed new state. Subtracting the exact linear equation gives the state error as minus the inverse of A applied to that residual. For a finite Hermitian Hamiltonian and the ordinary two-norm, every singular value of A is at least one, so its inverse norm is at most one. This gives a clean one-step absolute error bound for this specific system, without assuming that all unrelated implicit solves share the same conditioning.

Across many steps, solver residuals can accumulate and compete with truncation error. A tolerance fixed at an inconvenient absolute scale may dominate a timestep-refinement study. Renormalizing after every step can hide norm drift without removing phase error or replacing the missing solve accuracy. Adaptive tests should therefore compare step-subdivided states or selected observables and track both linear residuals and temporal changes. Small norm drift is useful evidence, but never the whole error estimate.

r=b−Ac^,c^−c=−A−1r,∥c^−c∥2≤∥r∥2(H=H†)r=b-A\widehat c,\qquad \widehat c-c=-A^{-1}r,\qquad \|\widehat c-c\|_2\le\|r\|_2\quad(H=H^\dagger)
The bound is for the stated Cayley linear system in an orthonormal finite basis. It controls one solve error, not the continuum spatial approximation or an arbitrary long-time phase error.↔ Scroll wide equations horizontally.

A norm-preserving but phase-biased two-level step

Take H=ϵdiag⁡(0,2)H=\epsilon\operatorname{diag}(0,2), initial plus state, and ϵΔt/ℏ=1\epsilon\Delta t/\hbar=1. Compute the Crank–Nicolson relative factor, final Bloch vector and excited-state probability. Compare the numerical relative phase with the exact one.

  1. The ground-state factor is one. For the excited level the ratio is (1−i)/(1+i)=−i(1-i)/(1+i)=-i, so the final state is (∣0⟩−i∣1⟩)/2(|0\rangle-i|1\rangle)/\sqrt2.

  2. Both probabilities remain one half, and the Bloch vector is (0,−1,0)(0,-1,0). The state norm is exactly one under this solved finite update.

  3. Its relative phase is minus pi over two, while exact evolution gives minus two radians. The phase discrepancy is 2−π/2≈0.42922-\pi/2\approx0.4292 radians in magnitude, even though norm, energy and these basis probabilities are conserved.

The result

The Cayley step produces a minus-i relative factor and Bloch y minus one, but its relative phase differs from exact evolution by about 0.4292 radians.

Explore the model

Interactive lab

Interactive lab · quantum

Equal populations can hide a phase difference

Reset and set mixing angle ninety degrees, phase minus ninety degrees; predict Bloch coordinates (0,−1,0)(0,-1,0) and equal basis probabilities. Compare a fresh Reset run at ninety degrees, phase zero; predict Bloch x one with the same probabilities. This pure-state geometry does not execute Crank–Nicolson or measure a solver residual.

Equal populations can hide a phase differenceA normalized two-level pure state, not an interacting field simulation. Playback adds one radian per second to the relative phase; basis probabilities stay fixed. Current readings: Probability |0⟩: 71.1%; Probability |1⟩: 28.9%; Bloch (x, y, z): 0.78, 0.45, 0.42. All parameters are adjustable with the labeled controls below.xyz|0⟩|1⟩Basis probabilitiesP(0)71.1%P(1)28.9%
Bloch vector / P(0)P(1)Phase latitude
Probability |0⟩
71.1%
Probability |1⟩
28.9%
Bloch (x, y, z)
0.78, 0.45, 0.42
t = 0 s

A normalized two-level pure state, not an interacting field simulation. Playback adds one radian per second to the relative phase; basis probabilities stay fixed. Parameter changes pause playback and recompute the model at the current elapsed time. Reset restores both parameters and returns time to zero. Playback starts paused. Adjust either parameter, then choose Play or Step.

Look forPopulation agreement alone does not control coherent phase error.

Check your understanding

For a solved Crank–Nicolson update with constant Hermitian H, which error can survive exact norm and energy conservation?

Try it yourself

For the scalar Cayley system A=1+iA=1+i, b=1−ib=1-i, suppose a solver returns c^=0.01−i\widehat c=0.01-i. Compute its residual, its error relative to the exact solution, and its squared norm. Does renormalization make the linear equation exact?

Take a moment to make a prediction or work it out on paper. Open a hint when you need one.

Give me a hint

The exact ratio is minus i. Compute b minus A times the supplied approximation before taking magnitudes.

Compare with the solution

The residual is r=−0.01(1+i)r=-0.01(1+i) with magnitude 0.0120.01\sqrt2. The exact solution is minus i, so the state error magnitude is 0.01, equal to residual magnitude divided by ∣A∣=2|A|=\sqrt2. The approximate squared norm is 1.00011.0001. Renormalizing can restore unit norm but changes the supplied vector without generally satisfying the original linear equation or removing its phase error. Solve accuracy must be checked through the residual and its propagation effect.

Take it with you

  • The Cayley proof needs finite Hermitian structure and an accurate solve.
  • Unit modulus and conserved energy do not guarantee the correct relative phase.
  • Linear residuals belong in the same budget as timestep and basis errors.

Read further

Original sources and deeper explanations to follow at your own pace.

Optional progress, saved in this browser when available. Revisit any lesson at any time.

Lesson not marked complete.