  The jump unravelling of a quantum master equation decomposes the dynamics of
an open quantum system into abrupt jumps, interspersed by periods of coherent
dynamics where no jumps occur. Simulating these jump trajectories is
computationally expensive, as it requires very small time steps to ensure
convergence. This computational challenge is aggravated in regimes where the
coherent, Hamiltonian dynamics are fast compared to the dissipative dynamics
responsible for the jumps. Here, we present a quantum version of the Gillespie
algorithm that bypasses this issue by directly constructing the waiting time
distribution for the next jump to occur. In effect, this avoids the need for
timestep discretisation altogether, instead evolving the system continuously
from one jump to the next. We describe the algorithm in detail and discuss
relevant limiting cases. To illustrate it we include four example applications
of increasing physical complexity. These additionally serve to compare the
performance of the algorithm to alternative approaches -- namely, the
widely-used routines contained in the powerful Python library QuTip. We find
significant gains in efficiency for our algorithm and discuss in which regimes
these are most pronounced. Publicly available implementations of our code are
provided in Julia and Mathematica.
