Large-scale stochastic simulation of open quantum systems

TL;DR

This work tackles the challenge of simulating large-scale open quantum systems governed by Lindblad dynamics. It introduces the tensor jump method (TJM), a scalable trajectory-based approach that blends Monte Carlo wave function unravelings with matrix product state techniques, incorporating a dynamic TDVP and a second-order Strang splitting to suppress time-step errors. A novel sampling MPS enables trajectory sampling at arbitrary times without inflating costs, and a site-local dissipative contraction keeps entanglement growth in check. The authors demonstrate strong convergence guarantees and practical speedups, achieving accurate simulations up to 1000 spins on a consumer CPU, with benchmarking against MPO solvers and analytical steady-state results. This method promises substantial utility for benchmarking, design, and exploration of noisy quantum dynamics in regimes previously beyond classical capabilities, potentially informing stable quantum hardware development and large-scale quantum simulations.

Abstract

Understanding the precise interaction mechanisms between quantum systems and their environment is crucial for advancing stable quantum technologies, designing reliable experimental frameworks, and building accurate models of real-world phenomena. However, simulating open quantum systems, which feature complex non-unitary dynamics, poses significant computational challenges that require innovative methods to overcome. In this work, we introduce the tensor jump method (TJM), a scalable, embarrassingly parallel algorithm for stochastically simulating large-scale open quantum systems, specifically Markovian dynamics captured by Lindbladians. This method is built on three core principles where, in particular, we extend the Monte Carlo wave function (MCWF) method to matrix product states, use a dynamic time-dependent variational principle (TDVP) to significantly reduce errors during time evolution, and introduce what we call a sampling MPS to drastically reduce the dependence on the simulation's time step size. We demonstrate that this method scales more effectively than previous methods and ensures convergence to the Lindbladian solution independent of system size, which we show both rigorously and numerically. Finally, we provide evidence of its utility by simulating Lindbladian dynamics of XXX Heisenberg models up to a thousand spins using a consumer-grade CPU. This work represents a significant step forward in the simulation of large-scale open quantum systems, with the potential to enable discoveries across various domains of quantum physics, particularly those where the environment plays a fundamental role, and to both dequantize and facilitate the development of more stable quantum hardware.

Figures (10)

• Figure 1: Algorithms of the Monte Carlo wave function (MCWF) method and the tensor jump method (TJM). (a) The MCWF algorithm stochastically simulates quantum trajectories by evolving a wavefunction $\ket{\Psi(t)}$ with a non-Hermitian effective Hamiltonian $e^{-i H^{\text{eff}} \delta t}$. A stochastic norm loss $\delta p = 1 - \langle \Psi^{(i)}(t+\delta t)|\Psi^{(i)}(t+\delta t)\rangle$ is compared against a uniform random number $\epsilon \in [0,1]$. If $\epsilon \geq \delta p$, the result of the process is the dissipated state $|\Psi^{(i)}(t+\delta t)\rangle$. Otherwise, a quantum jump operator $L_j$ is applied based on computed jump probabilities to the pre-evolved state. The output of this stochastic process is then normalized. This process is repeated iteratively to simulate dissipation until the final time $T$. (b) The TJM extends this idea to tensor networks. A sampling MPS $\ket{\Phi(j\delta t)}$ is evolved by alternating applications of the dissipative sweep $\mathcal{D}$, a stochastic jump process $\mathcal{J}_\epsilon$, and dynamic TDVP $\mathcal{U}$. At any point, the quantum trajectory $\ket{\Psi(j \delta t)}$ can be sampled from $\Phi$ by applying a final sampling layer. The method retains the structure of MCWF while enabling low timestep error from second-order Trotterization without additional overhead. Further details of the operators $\mathcal{U}$, $\mathcal{D}$, and $\mathcal{J}_\epsilon$ are provided in \ref{['sec:TJM']}.
• Figure 2: This figure visualizes each step of the various forward- and backward-evolving equations of 1TDVP and 2TDVP indicated by the black dashed line. This shows the reduced practical form of the network caused by the MPS's mixed canonical form. The tensors surrounded by the dashed lines (corresponding to $H^{\text{eff}}$) are contracted, exponentiated with the Lanczos method, then applied to the remaining tensors to update them. (Top) 1TDVP forward-evolving and 2TDVP backward-evolving network where $H^\text{eff}_3$ ($\tilde{H}^\text{eff}_{3,4}$) is a degree-6 tensor. (Middle) 1TDVP backward-evolving network where $\tilde{H}^\text{eff}_3$ is a degree-4 tensor. (Bottom) 2TDVP forward-evolving network where $H^\text{eff}_{3,4}$ is a degree-8 tensor.
• Figure 3: Two components of the tensor jump method (TJM). (Top) Illustration of the application of the factorized dissipative MPO $\mathcal{D}$ (constructed from local matrices) to an MPS $\ket{\Psi}$. Each local tensor corresponds to the exponentiation of local jump operators: $\mathcal{D}_\ell = e^{-\frac{1}{2} \delta t \sum_{j \in S(\ell)} L_j^{[\ell] \dagger} L_j^{[\ell]}}$. This operation does not change the bond dimension and does not require canonicalization. (Bottom) Visualization of the tensor network required to compute $\delta p_m$ for a given jump operator $L_m$, which corresponds to the expectation value $\langle L^\dagger_m L_m \rangle$. By putting the MPS in mixed canonical form centered at site $\ell$, contractions over tensors to the left and right reduce to identities. This enables efficient computation of the probability distribution $\Pi(t)$ via a sweep across the network, evaluating local jump probabilities $L_j$ at each site $j \in S(\ell)$.
• Figure 4: Convergence benchmarks of the TJM. (a) Error in the local observable $\langle X^{[5]}\rangle$ at time $Jt = 1$ as a function of the number of trajectories $N$, for several time step sizes $\delta t$. Each point represents the average over 1000 batches. Dotted lines show the corresponding first-order Trotterization results for comparison. The shaded region indicates the standard deviation for $\delta t = 0.1$. The error follows the predicted scaling $\sim C / \sqrt{N}$, with $C \approx 0.1$, demonstrating that the second-order Trotter scheme yields low time-step errors such that $N$ dominates over $\delta t$ in determining convergence. (b) Convergence of the two-site correlator $\langle X^{[i]} X^{[i+1]}\rangle$ over time $Jt \in [0, 10]$ at fixed $\delta t = 0.1$, shown as a function of trajectory number $N$ and bond dimension $\chi$. Each panel displays the local observable error, averaged over 1000 batches of $N$ trajectories. The colormap is centered at an error threshold $\epsilon = 10^{-2}$: blue indicates lower error, red higher. While increasing $N$ significantly improves global accuracy, a larger bond dimension is still required to resolve localized dynamical features.
• Figure 5: Convergence of the TJM with increasing bond dimension $\chi = \{2, 4, 8\}$ for a 30-site noisy XXX Heisenberg model with parameters $J = h = 1$ and a domain-wall initial state (wall at site 15). The noise model includes local relaxation and excitation channels $\sigma_\pm$ with equal strength $\gamma_\pm = 0.1$. (a–d) Site-resolved magnetization $\langle Z^{[i]} \rangle$ over time, compared against a reference simulation using an MPO with bond dimension capped at $D = 400$. (e–h) Time evolution of $\langle Z^{[i]} \rangle$ for selected sites $i = 1, 5, 10, 15$, showing quantitative agreement between TJM and MPO as $\chi$ increases. All simulations use a timestep $\delta t = 0.1$, and TJM is averaged over $N = 100$ trajectories. While the MPO simulation requires over 24 hours, the TJM runs complete in under 5 minutes. This highlights the ability of TJM to efficiently and accurately reproduce both global and local dynamics at modest bond dimension, making it a practical tool for rapid prototyping of noisy quantum dynamics.

Theorems & Definitions (11)

• Definition 1: Density matrix variance
• Theorem 2: Convergence of TJM
• proof
• Theorem 3: Equivalence of MCWF and Lindblad master equation
• proof
• Lemma 4: Frobenius variance
• Definition 5: Density matrix variance and standard deviation
• Lemma 6: Frobenius variance
• proof
• Theorem 7: Convergence of TJM