The simulation of dissipative quantum dynamics [1]} is a problem which arises in diverse areas of physics, including ultrafast spectroscopy, chemical physics, quantum optics, quantum biology, quantum computing, and quantum information technology [2]}, [3]}, [4]}, [5]}, [6]}. Some problems involving open quantum systems in the presence of a driving field, which include optimal control of open quantum systems [7]}, [8]} and simulating coherent pulse propagation using the Maxwell-Schrödinger equations [9]} when accounting for divergence- and propagation-induced decoherence [10]}, [11]}, [12]}, require simulating the system iteratively for different waveforms of the field. Many methods for simulating open quantum systems are available, including Monte Carlo methods [13]}, time evolving density matrix methods, and path integral approaches. However, they often become unsuitable for iterative methods by making any computations involving them undesirably expensive, and, for optimal control problems, by not having a straightforward way of computing gradients.
