The evolution of open quantum systems is a fundamental topic in various scientific fields. During time propagation, the environment occasionally makes measurements, forcing the system's wave function to collapse randomly. The von Neumann density matrix incorporates the statistics involved in these random processes, and its time development is often described by Markovian quantum master equations that incorporate a dissipator. For large systems, the complexity of the dissipator grows with the increasing number of possible measurements, posing conceptual and severe computational challenges. This article introduces a stochastic representation of the dissipator, using bundled measurement operators to address this complexity. Taking the Morse oscillator as an example, we demonstrate that small samples of bundled operators capture the system's dynamics. This stochastic bundling is different from the stochastic unraveling and the jump operator formalism and offers a new way of understanding quantum dissipation and decoherence.

Abstract Simulating mixed-state evolution in open quantum systems is crucial for various chemical physics, quantum optics, and computer science applications. These simulations typically follow the Lindblad master equation dynamics. An alternative approach known as quantum state diffusion unraveling is based on the trajectories of pure states generated by random wave functions, which evolve according to a nonlinear Itô-Schrödinger equation (ISE). This study introduces weak first- and second-order solvers for the ISE based on directly applying the Itô-Taylor expansion with exact derivatives in the interaction picture. We tested the method on free and driven Morse oscillators coupled to a thermal environment and found that both orders allowed practical estimation with a few dozen iterations. The variance was relatively small compared to the linear unraveling and did not grow with time. The second-order solver delivers much higher accuracy and stability with bigger time steps than the first-order scheme, with a small additional workload. However, the second-order algorithm has quadratic complexity with the number of Lindblad operators as opposed to the linear complexity of the first-order algorithm.