Gillespie algorithm for quantum jump trajectories

The jump unraveling of a quantum master equation decomposes the dynamics of an open quantum system into abrupt jumps, interspersed by periods of coherent dynamics when no jumps occur. Such open quantum systems are ubiquitous in quantum optics and mesoscopic physics, hence the need for efficient techniques for their stochastic simulation. Numerical simulation techniques fall into two main categories. The first splits the evolution into small time steps and determines stochastically for each step if a jump occurs or not. The second, known as Monte Carlo wave-function simulation, is based on the reduction of the norm of an initially pure state in the conditional no-jump evolution. It exploits the fact that the purity of the state is preserved by the finest unraveling of the master equation. In this paper, we present an alternative method for the simulation of the quantum jump unraveling, inspired by the classical Gillespie algorithm. The method is particularly well suited for situations in which a large number of trajectories is required for relatively small systems. It allows for non-purity-preserving dynamics, such as the ones generated by partial monitoring and channel merging. We describe the algorithm in detail and discuss relevant limiting cases. To illustrate it, we include four example applications of increasing physical complexity and discuss the performance of the algorithm across regimes of interest for open quantum systems simulation. Publicly available implementations of our code are provided in julia and mathematica.