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 $\Delta 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}(\Delta t)$ and $\mathcal{O}(\Delta t^2)$ 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^3 2^{2K_{\max}} K_{\max}^2)$ and a space complexity of $\mathcal{O}(M^2 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 ($M\geqslant 3$). Numerical experiments are carried out to illustrate the validity and efficiency of our method.