# Gaussian Quantum Trajectories for the Variational Simulation of Open Quantum-Optical Systems
## 1 Introduction
We are interested in a method that is variational and provides a description on the level of trajectories. Some current variational (Time-Dependent Variational Principle (TDVP) ([PERSON], 1985)) approaches to the description of open systems aim to describe the system at the level of the master equation ([PERSON], 1985; 1988; 1988; 1988), with ansatzes including the Gutzwiller density matrix-([PERSON], 1990; 1991; 1992), matrix-product state ([PERSON], 1990; 1992) and matrix-product operator ([PERSON], 1990; 1992; 1993; 1994) methods. Variational descriptions on the trajectory level have so far focused on lattice-like systems: systems where the complexity arises mainly due to a large amount of modes, while the Hilbert space per mode remains small, with methods such as Gutzwiller Monte Carlo ([PERSON], 1990; 1992) or based on t-DMRG ([PERSON], 1992; 1993; 1994). Somewhat surprisingly, except for some works on the back-reaction of measurements on quantum many-body systems ([PERSON], 1990; 1992), the class of Gaussian states ([PERSON], 1990) has not received much attention as a variational ansatz for the simulation of quantum trajectories. This stands in stark contrast with equilibrium systems, where the Gaussian ansatz has been fruitfully exploited, regarding Hartree-Fock-Bogoliubov ([PERSON], 1963; [PERSON], 1963), and with the important application of time-dependent BCS/BdG ([PERSON], 1963; 1964) theory. It is the purpose of this work to develop the formalism of a Gaussian variational description of quantum trajectories. Since the Wigner distribution of a Gaussian state is entirely positive, this description is expected to work well in the regime of many particles and weak interactions, when the Wigner distribution is positive everywhere ([PERSON], 1992). In this regime, a different stochastic technique, the so-called Truncated Wigner Approximation (TWA), is widely used ([PERSON], 1992; 1993; 1994; 1995; 1996; 1997; 1998; 1999;  
'Gaussian states', to which we will refer for clarity as '\(XP\)-Gaussian states'. We will apply such an \(XP\)-Gaussian trajectory method to study the stationary state of a driven-dissipative Kerr cavity in the bistability regime [59]. Typical Bogoliubov expansion around the mean-field value is fundamentally unable to describe the hopping between the two branches. We will see that these \(XP\)-Gaussian methods on the other hand can make reasonable predictions of the exact solution. In Section 4, however, we will see that \(XP\)-Gaussian methods are unable to describe situations where a large variance in phase is present. A convenient alternative that does succeed in providing accurate descriptions is comprised of the \(N\Theta\)-Gaussian methods, meaning the assumption that the state of the trajectory is Gaussian in density and phase. In Section 5, we will address computational issues and compare the performance with TWA and exact trajectory methods. Finally, in Section 6, we conclude the work. Throughout our work, we will investigate the effect of different unraveling schemes (photon counting and heterodyne detection) on the accuracy of our variational method. Code S1 used for the numerical simulations is available as Supplementary Material.