Classical simulation of non-Gaussian bosonic circuits

We propose efficient classical algorithms which (strongly) simulate the action of bosonic linear optics circuits applied to superpositions of Gaussian states. Our approach relies on an augmented covariance matrix formalism to keep track of relative phases between individual terms in a linear combination. This yields an exact simulation algorithm whose runtime is polynomial in the number of modes and the size of the circuit and quadratic in the number of terms in the superposition. We also present a faster approximate randomized algorithm whose runtime is linear in this number. Our main building blocks are a formula for the triple overlap of three Gaussian states and a fast algorithm for estimating the norm of a superposition of Gaussian states up to a multiplicative error. Our construction borrows from earlier work on simulating quantum circuits in finite-dimensional settings, including, in particular, fermionic linear optics with non-Gaussian initial states and Clifford computations with nonstabilizer initial states. It provides algorithmic access to a practically relevant family of non-Gaussian bosonic circuits.