An efficient algorithm to compute the exponential of skew-Hermitian matrices for the time integration of the Schr\"odinger equation

TL;DR

This paper introduces a computationally efficient algorithm for approximating the exponential of skew-Hermitian matrices, crucial for quantum mechanics simulations, using Chebyshev polynomials with reduced matrix multiplications.

Contribution

It develops a novel Chebyshev polynomial-based method with optimized degrees and matrix operations, outperforming traditional rational Padé and Taylor schemes in efficiency.

Findings

More efficient than Padé and Taylor methods across tolerances

Uses fewer matrix-matrix products for high-degree polynomials

Suitable for exponential integrators in Schrödinger equation simulations

Abstract

We present a practical algorithm to approximate the exponential of skew-Hermitian matrices up to round-off error based on an efficient computation of Chebyshev polynomials of matrices and the corresponding error analysis. It is based on Chebyshev polynomials of degrees 2, 4, 8, 12 and 18 which are computed with only 1, 2, 3, 4 and 5 matrix-matrix products, respectively. For problems of the form $\exp(-iA)$, with $A$ a real and symmetric matrix, an improved version is presented that computes the sine and cosine of $A$ with a reduced computational cost. The theoretical analysis, supported by numerical experiments, indicates that the new methods are more efficient than schemes based on rational Pad\'e approximants and Taylor polynomials for all tolerances and time interval lengths. The new procedure is particularly recommended to be used in conjunction with exponential integrators for the numerical time integration of the Schr\"odinger equation.