Second-order discretization of Dyson series: iterative method, numerical analysis and applications in open quantum systems

We propose a general strategy to discretize the Dyson series without applying direct numerical quadrature to high-dimensional integrals, and extend this framework to open quantum systems. The resulting discretization can also be interpreted as a Strang splitting combined with a Taylor expansion. Based on this formulation, we develop a numerically exact iterative method for simulation system-bath dynamics. We propose two numerical schemes, which are first-order and second-order in time step Δt respectively. We perform a rigorous numerical analysis to establish the convergence orders of both schemes, proving that the global error decreases as mathcal{O}(Δt) and mathcal{O}\(Δt²\) for the first- and second-order methods, respectively. In the second-order scheme, we can safely omitted most terms arising from the Strang splitting and Taylor expansion while maintaining second-order accuracy, leading to a substantial reduction in computational complexity. For the second-order method, we achieves a time complexity of mathcal{O}\(M³ 2^2K_max K_max²\) and a space complexity of mathcal{O}\(M² 2^2K_max K_max\) where M denotes the number of system levels and K_max the number of time steps within the memory length. Compared with existing methods, our approach requires substantially less memory and computational effort for multilevel systems $Mgeqslant 3$. Numerical experiments are carried out to illustrate the validity and efficiency of our method.