A Posteriori Error Analysis and Rank-Adaptive Structure-Preserving Numerical Integration of Large-Scale Lindblad Matrix Differential Equations

Abstract

Large-scale Lindblad matrix differential equations are difficult to integrate because a dense density matrix requires quadratic storage while standard time integrators can destroy positivity. This paper develops a second-order rank-adaptive structure-preserving integrating-factor method, denoted RA-SPIF2, together with a fully computable a posteriori error analysis in the trace norm. The update is derived from an explicit midpoint rule in the interaction picture and is implemented as a sum of positive semidefinite congruences acting directly on a tall factor. A thin QR factorization and a small Hermitian eigenproblem select the numerical rank. The discarded eigenvalue mass yields an exact trace-distance identity for the compression perturbation, while normalization enforces unit trace. For error control, two accepted half steps define a quadratic continuous reconstruction. Contractivity of the exact Lindblad semigroup gives a global residual bound, and a Bernstein representation converts each local residual integral into a certified expression involving only three trace norms. These norms are evaluated through small projected Hermitian matrices, without forming the density matrix or the Liouvillian. Numerical experiments for a driven thermal oscillator exhibit observed orders 1.74, 1.87, and 1.93 under refinement; the certified estimator remains reliable with effectivities from 22.36 to 17.49, while a Gauss residual estimator gives affectivities from 7.20 to 5.51. Rank studies reveal a clear temporal-error plateau, and coarse-step comparisons show positivity loss for classical RK4 but only roundoff-level violations for RA-SPIF2. A boundary-driven spin chain is advanced to dimension 16,384 with rank 10, using 2.5 MiB for the state factor instead of 4,096 MiB for a dense complex matrix.