Precision Diagnostics in Quantum Calculations

Numerical precision is adequate only when the observable of interest is stable at the required scale. Energy eigenvalues may each have many accurate digits while their splitting is unresolved; a small residual may reflect an accurate solve of an inadequate grid. Controlled changes of precision, discretization, and solver tolerance help distinguish these cases.

Helpful background. Floating-point arithmetic explains representation, unit roundoff, and cancellation. Numerical units explains scaling; array representations fixes the inner products being computed.

Energy scales and arithmetic range

Use $u$ for the unit roundoff, equal to $2^{-53}\approx1.11\times10^{-16}$ for binary64 rounding to nearest. The local rounding model bounds individual basic operations under its stated conditions; it does not bound an entire calculation’s final error.

Absolute spacing grows with magnitude, while relative spacing remains approximately constant for normal numbers. Near $2^{53}$, not every consecutive integer is representable in binary64. A small energy correction added to a sufficiently large reference energy can disappear. Shifting a Hamiltonian by a scalar multiple of the identity may improve the scale of intermediate expressions, but computing the shift from already inaccurate matrix elements cannot recover lost information.

Subnormal numbers extend the range near zero with reduced relative precision. Overflow produces infinities, while invalid operations can produce NaNs. Check these explicitly: a convergence test involving NaN may evaluate to false and accidentally terminate a loop. Hardware modes that flush subnormals to zero can alter very small amplitudes.

Exponentials require particular care. A tunneling weight $P=Ae^{-B}$ can underflow even when $\log P=\log A-B$ is accurately representable. Report and test the exponent and prefactor separately; a printed zero does not establish a vanishing physical probability. For thermal weights, subtract a common reference energy before exponentiating and retain the removed factor when computing an absolute partition function. Thermal-state numerics develops this normalization.

The magnitude of $e^{-iEt/\hbar}$ does not underflow, but forming an enormous phase argument can lose the low-order information that determines its value modulo $2\pi$. Factoring out a known common phase or reference energy can help before the information is lost. Rescaling a wavefunction can likewise prevent overflow or underflow, provided the scale factors needed for the requested observable are retained.

Cancellation in level splittings

Suppose $E_1$ and $E_2$ are each of order $E_*$ and known with absolute error of order $uE_*$. Their difference $\Delta=E_2-E_1$ can have relative error of order

$$\frac{uE_*}{|\Delta|}.$$

If $E_*=1$ and $\Delta=10^{-12}$, this scale is about $10^{-4}$ even though each energy has many accurate digits. The subtraction may itself be exactly rounded; the difficulty is the uncertainty already present in its operands. Near-degenerate tunneling splittings therefore deserve a separate error budget, symmetry-aware formulation, and precision study.

Algebraic reformulation helps when the small quantity can be computed directly, as the rationalized square-root example illustrates. For a normalized state close to an energy eigenstate, use the Rayleigh-quotient residual to evaluate its energy variance as $\|(H-\langle H\rangle)\psi\|^2$, instead of subtracting nearly equal values of $\langle H^2\rangle$ and $\langle H\rangle^2$. This avoids one cancellation mechanism; the matrix action and energy expectation still need adequate accuracy.

Compare a claimed splitting with the original eigenpair residuals, changes under basis or grid refinement, and the symmetry that predicts or forbids the splitting. For clustered levels, compare spectral projectors as well as energies: numerical eigenvector rotation inside a nearly degenerate subspace need not represent a physical change of state.

Summation, normalization, and orthogonality

Sequentially summing $n$ terms has a standard absolute error bound involving

$$\gamma_{n-1}\sum_j|a_j|, \qquad \gamma_k=\frac{ku}{1-ku},\quad ku<1.$$

The relative error can be large when the exact sum is small. Inner products of nearly orthogonal states have precisely this cancellation pattern. Pairwise or compensated summation can reduce accumulation error, but cannot remove uncertainty in the input coefficients. Complex inner products require conjugation independently of summation accuracy.

A naive norm $\sqrt{\sum_j|c_j|^2}$ can overflow even when the final norm is representable. Scale first by $s=\max_j|c_j|$, compute $s\sqrt{\sum_j|c_j/s|^2}$, and treat $s=0$ separately. Established linear-algebra libraries use related safe scaling strategies.

Loss of orthogonality in a long Krylov recurrence can create repeated apparent eigenvalues. Checking only energy residuals misses this issue; measure the Gram matrix of computed vectors as well. Reorthogonalization addresses an algorithmic accumulation mechanism, whereas increasing precision merely changes its scale.

Refinement can amplify roundoff

For a smooth function on a dimensionless coordinate grid, the centered second derivative is

$$D_h^{(2)}f(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}.$$

Taylor expansion gives truncation error $h^2f^{(4)}(x)/12+O(h^4)$. If evaluations of the numerator’s terms have absolute errors of order $u$ on a fixed function scale, division by $h^2$ amplifies those errors. A useful schematic envelope is

$$\mathcal E(h)\approx C_1h^2+C_2\frac{u}{h^2}, \qquad C_1,C_2>0.$$

The constants depend on the function, evaluation method, and local scale; this is not an exact error formula for every grid eigenvalue. For a dimensional coordinate, the constants carry the corresponding units. Under the stated model the best balance satisfies

$$h_*^4=\frac{C_2u}{C_1}, \qquad \mathcal E(h_*)=2\sqrt{C_1C_2u}.$$

Refining past this region can make the result worse. Increasing precision shifts the balance toward smaller $h$, whereas enlarging a finite box addresses a separate boundary error. The finite-difference Schrödinger calculation separates these limits from algebraic eigensolver accuracy.

Numerical repairs and physical constraints

A small Hermiticity defect can arise from rounding in assembly. Projecting $H$ to $(H+H^\dagger)/2$ is defensible only when Hermiticity is an exact property of the intended model and the correction is consistent with the assembly-error budget. Measure the defect first and retain that diagnostic. Complex magnetic-field matrix elements need not violate Hermiticity; a deliberate absorbing potential does, and symmetrizing it changes the model.

The same principle applies to normalizing a state or projecting a density operator onto positive matrices. A repaired invariant cannot independently validate the step that was repaired. Record the defect before correction and test observables against an independent benchmark or refinement study.

A precision experiment

Hold the model, box, grid, and solver criteria fixed while increasing arithmetic precision. Compare a target observable, residual, and invariant. Then vary the discretization separately. If increasing precision changes nothing but grid refinement does, discretization dominates. If finer grids stop improving and higher precision restores the expected trend, roundoff is a plausible limiting mechanism.

A reproducible calculation need not be bitwise identical across parallel reduction orders. Record the target observable, units, grid or basis, precision, solver and library, tolerances, and relevant reduction policy. Compare residuals and benchmark errors with scientifically meaningful tolerances. A dramatic discrepancy is not explained merely by saying that floating-point addition is nonassociative.

Exercises

A pair of computed energies is accurate to an absolute tolerance $10^{-9}E_*$ each. Give a conservative absolute uncertainty bound for their difference and assess a claimed splitting of $10^{-10}E_*$.

Solution

The triangle inequality gives $2\times10^{-9}E_*$. The claimed splitting is smaller than this bound and is unresolved by those energy estimates. Correlated errors might cancel, but their correlation must be established independently before using a smaller uncertainty.

Why can renormalizing a wavefunction after every time step hide a faulty propagator?

Solution

Renormalization forces the norm diagnostic to pass while leaving phase, population, and relative-amplitude errors. Record the norm before correction and test the unmodified propagation against an analytic solution or step refinement.

Use the residual identity to explain why directly subtracting two energy moments is a poor variance diagnostic near an eigenstate. Does a zero residual identify the ground state?

Solution

The moments $\langle H^2\rangle$ and $\langle H\rangle^2$ approach the same value, so uncertainty in either can dominate their small difference. The equivalent squared residual norm computes a small vector before taking its norm and avoids that final subtraction. Every exact energy eigenstate has zero residual; the diagnostic alone does not identify the ground state. The matrix-diagonalization exercise proves the identity.

Minimize the schematic error $\mathcal E(h)=C_1h^2+C_2u/h^2$. If $u$ is reduced by a factor of $16$, how do the optimal spacing and minimum envelope change?

Solution

Setting $\mathcal E'(h)=2C_1h-2C_2u/h^3=0$ gives $h_*=(C_2u/C_1)^{1/4}$. At this point both error terms equal $\sqrt{C_1C_2u}$. Reducing $u$ by $16$ halves $h_*$ and reduces the minimum envelope by $4$. This prediction assumes the constants and evaluation procedure are otherwise fixed; an actual error curve can oscillate or have additional contributions.

A matrix has a small anti-Hermitian part. What must be established before replacing it by its Hermitian part?

Solution

Establish that the intended Hamiltonian is Hermitian, that the defect is explained by the numerical assembly rather than a missing conjugate or boundary term, and that its size is compatible with the desired observable accuracy. Report the correction. If the anti-Hermitian part represents absorption or decay, removing it changes the physical problem. A defect small relative to the overall matrix norm can still matter for a much smaller spectral splitting.

References

• Goldberg, D., “What Every Computer Scientist Should Know About Floating-Point Arithmetic,” ACM Computing Surveys 23, 5–48 (1991), authorized reprint — rounding, cancellation, exceptional values, and arithmetic semantics.
• Higham, N. J., Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, publisher edition — summation bounds and stable norm computation.
• Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007 — balancing truncation and roundoff in numerical differentiation.