Generic Method for Integrating Lindblad Master Equations

TL;DR

This paper introduces a generic Taylor-series-based method for integrating Lindblad master equations without vectorizing the density matrix, yielding significant memory and time savings for large open quantum systems and enabling seamless tensor-network implementations. By expanding the Lindbladian exponential and applying it directly to the initial state, the method achieves ${\cal O}(d^2)$ space and ${\cal O}(d^3)$ time complexity, with controllable truncation error. Through two illustrative examples—a driven two-level system and a dissipative Heisenberg chain—the authors demonstrate numerical accuracy, trace preservation, and compatibility with MPDO/MPO representations and METTS-based thermal states. Performance benchmarks show clear advantages over the standard vectorization approach and competitive performance relative to RK4, especially for large systems. The work suggests future extensions using Krylov subspaces and Magnus expansions to further enhance efficiency for time-dependent open-system dynamics, positioning the Taylor-series method as a practical, scalable solver for nonequilibrium quantum dynamics.

Abstract

The time evolution of Markovian open quantum systems is governed by Lindblad master equations, whose solution can be formally written as the Lindbladian exponential acting on the initial density matrix. By expanding this Lindbladian exponential into the Taylor series, we propose a generic method for integrating Lindblad master equations. In this method, the series is truncated, retaining a finite number of terms, and the iterative actions of Lindbladian on the density matrix follow the corresponding master equation. Our method offers significant improvements in numerical efficiencies both in memory cost and computation time, especially for systems with many degrees of freedom. Moreover, our proposed method can be integrated seamlessly with tensor networks. Two illustrative examples, a two-level system exhibiting damped Rabi oscillations and a driven dissipative Heisenberg chain, are used to demonstrate the validity of our method. The superiority of our method is benchmarked with detailed performance tests.

Figures (8)

• Figure 1: The decaying behavior of the relative error with the number of retained terms in the Taylor series method. The parameter $\Delta$ is given in terms of the Lindbladian norm and the time: $\Delta=||{\cal L}||t$.
• Figure 2: The time evolution of the density matrix of the dissipative two-level system initially in the excited state $\ket{1}$. The two diagonal elements are shown here. The solid lines are results calculated according to the vectorization method #1, the starred points according to TaylorSeries10, and the filled circles according to RK4. The parameter values adopted in numerical simulations are $E=\Omega=\hbar=1.0$ and $\Gamma=0.5$. In both panels, the time step $\delta t=0.1$ is taken for the vectorization method #1. For the other two methods, the time step $\delta t=0.5$ is taken in the left panel, and $\delta t=1.0$ in the right panel.
• Figure 3: Quantum trajectories of the matrix product state (MPS). The ensemble-averaged evolution of the system's density matrix can be unraveled in terms of individual quantum trajectories, each consisting of abrupt jumps occurring at random times and non-Hermitian evolution in between them. The effective non-Hermitian Hamiltonian is defined by $H_{\rm eff}=H-\frac{{\rm i}\hbar}{2}\sum_iL_i^{\dagger}L_i$. The effective non-Hermitian evolution operator ${\rm e}^{-{\rm i}H_{\rm eff}\delta t/\hbar}$ is represented as a matrix product operator (MPO; circumscribed by the dashed line) that is constructed from truncated Taylor series, in the same spirit as the Taylor series method. The quantum jump as sketched in the figure corresponds to a spin lowering operator acting on the rightmost site. The density matrix is represented by matrix product density operator (MPDO, or simply MPO), which can be constructed from individual quantum trajectories.
• Figure 4: The time evolution of the density matrix of the dissipative Heisenberg chain, which consists of $5$ spins. The system starts its evolution from a thermal state $\rho_{\rm eq}(0)={\rm e}^{-\beta H}/{\rm Tr}[{\rm e}^{-\beta H}]$. Two diagonal elements $\rho_{1,1}\equiv\braket{\uparrow\uparrow\uparrow\downarrow\downarrow|\rho|\uparrow\uparrow\uparrow\downarrow\downarrow}$ and $\rho_{2,2}\equiv\braket{\uparrow\downarrow\downarrow\downarrow\downarrow|\rho|\uparrow\downarrow\downarrow\downarrow\downarrow}$ are shown as representives. The solid lines are results calculated from quantum trajectories, while the starred points are according to the TaylorSeries10. The parameter values $J=\Gamma=\hbar=\beta=1.0$ are adopted in numerical simulations. In the quantum jump method, the time step $\delta t=0.1$ is taken, and $N=1000$ stochastic trajectories are generated to give the average evolution behavior. In TaylorSeries10, the time step $\delta t=0.2$ is taken. The diagonal elements are calculated by first representing pure states (ket and bra) as matrix product states (MPSs), and then contracting with MPO representation of the density matrix from both the bottom and the top sides. In the tensor-network calculations, the default parameter value cutoff=1.0e-14 is used when an MPO is applied to an MPO or MPS as the singular value decompositions (SVDs) are performed. This default parameter is defined as a real number $\epsilon$ by $\sum_{n\in{\rm discarded}}\sigma_n^2/\sum_n\sigma_n^2<\epsilon$, where $\{\sigma_n\}$ are singular values. The minor disagreement between the results from the two methods is due to the fluctuations in stochastic quantum trajectories.
• Figure 5: Benchmarks measuring the time elapsed when performing time evolution of the dissipative Heisenberg chain in one step according to the Lindblad master equation. All operators are represented as matrices. The computer program is coded in julia in single-thread mode. The computing platform is a PC with an Intel Core i7-12700H processor, Ubuntu 24.04.1 LTS OS, and 16 GB RAM. In the left panel, comparisons are made between TaylorSeries10 and the vectorization method, both #1 and #2. For vectorization method #1, the largest number of sites of the system that can be handled by the computer program is $L=6$. When $L=7$, the computer program was killed after long-time running, and when $L\ge 8$, the out-of-memory error occurred due to insufficient memory to store the Lindbladian ${\cal L}$ in the matrix representation. For vectorization method #2, the out-of-memory error also occurred due to the same reason when $L\ge 8$. In vectorization method #1, the built-in function exp is used for matrix exponential. As in TaylorSeries10, terms up to the 10-th order are retained in vectorization method #2. In the right panel, comparisons between TaylorSeries4 and RK4 are made and, as shown in the legend, the latter is implemented in different ways. Two optimization techniques are identified: one adopting the effective non-Hermitian Hamiltonian and the other using sparse matrices. In optimization scheme #1, only the first technique is applied. In optimization scheme #2, both techniques are applied. In both panels, the scaling behavior of computation time with the number of sites is observed, $\log_{10}({\rm time})\sim L\log_{10}8$, as expected.