Abstract
We present a new paradigm for the dynamical simulation of interacting many-boson open quantum systems. The method relies on a variational ansatz for the n-boson density matrix, in terms of a superposition of photon-added coherent states. It is most efficient for the simulation of driven-dissipative systems where the state is well described by quantum fluctuations on top of a displaced field, making it suitable for the simulation of several coupled modes with large occupation numbers, that are otherwise very challenging using a Fock-space expansion. We test our method on several examples, demonstrating its potential application to the predictive simulation of interacting bosonic systems and cat qubits.