# Gaussian Quantum Trajectories for the Variational Simulation of Open Quantum-Optical Systems
## 4 Phase Diffusion and the \(N\Theta\)-Gaussian Method
### N\(\Theta\)-Gaussian States  
The shape of the Wigner distribution suggests that a Gaussian description in terms of density and phase may provide a better approximation to the quantum trajectory wave functions. It must be noted that strictly speaking, a well-defined Hermitian phase operator does not exist [65]. Nevertheless, we will work with the Dirac definition \(\hat{\mathsf{a}}=:e^{i\hat{\theta}}\sqrt{\hat{\mathsf{n}}}\), where \(\hat{\theta}\) defines the phase and \(\hat{\mathsf{n}}\) is the familiar particle number operator. The practical use of this phase is paramount in quantum hydrodynamics [66]. This realm of use coincides with our approximation of density and phase to be continuous and unbounded operators, justified when the density is sufficiently high and the state is localized in phase-space, where the second demand allows avoiding complications arising from the multivaluedness of phase. Density and phase are conjugate variables: \(\left[\hat{\mathsf{n}},\hat{\theta}\right]=i\). Assuming the state to be Gaussian in \(\hat{\mathsf{n}}\) and \(\hat{\theta}\), the independent (real) expectation values to take into account are \(\left\langle\hat{\mathsf{n}}\right\rangle\), \(\left\langle\hat{\delta}_{n}\hat{\delta}_{n}\right\rangle=\left\langle\hat{ \mathsf{n}}\,\hat{\mathsf{n}}\right\rangle-\left\langle\hat{\mathsf{n}}\right\rangle ^{2}\), \(\left\langle\hat{\theta}\right\rangle\), \(\left\langle\hat{\delta}_{\theta}\hat{\delta}_{\theta}\right\rangle=\left\langle \hat{\theta}\hat{\theta}\right\rangle-\left\langle\hat{\theta}\right\rangle^{2}\) and \(\left\langle\hat{\delta}_{n}\hat{\delta}_{\theta}\right\rangle_{\text{sym}}= \left\langle\hat{\mathsf{n}}\,\hat{\theta}\right\rangle/2+\left\langle\hat{ \theta}\,\hat{\mathsf{n}}\right\rangle/2-\left\langle\hat{\mathsf{n}}\right\rangle \left\langle\hat{\theta}\right\rangle\). From these correlation functions, expectation values of products of annihilation and creation operators can be approximately computed by expanding the exponential and using [PERSON]'s theorem. For example,  
\[\left\langle\sqrt{\hat{\mathsf{n}}}e^{-i\hat{\theta}}\right\rangle\approx \sqrt{\left\langle\hat{\mathsf{n}}\right\rangle}e^{-i\left\langle \hat{\theta}\hat{\delta}_{n}\right\rangle} \tag{22}\] \[\times\left(1+\frac{1}{2\left\langle\hat{\mathsf{n}}\right\rangle }(\frac{1}{2}-i\left\langle\hat{\delta}_{n}\hat{\delta}_{\theta}\right\rangle_ {\text{sym}})-\frac{\left\langle\hat{\delta}_{n}\hat{\delta}_{n}\right\rangle }{8\left\langle\hat{\mathsf{n}}\right\rangle^{2}}\right).\]  
Here, phase-phase correlations were kept up to all orders, but phase-density correlations were truncated at second order. This approximation is again valid when the average particle number is sufficiently large. Using the given set of expectation values, derivation of the stochastic equations of motion, similarly to the case of \(XP\)-Gaussian states, proceeds as described in Section 2. For example, for the deterministic evolution under photon-counting, we obtain for the density:  
\[\partial_{t}\left\langle\hat{\mathsf{n}}\right\rangle=-\gamma\left\langle\hat {\delta}_{n}\hat{\delta}_{n}\right\rangle+2\operatorname{Im}\left[F\left\langle \sqrt{\hat{\mathsf{n}}}e^{-i\hat{\theta}}\right\rangle\right] \tag{23}\]  
while jumps:  
\[\left\langle\hat{\mathsf{n}}\right\rangle\xrightarrow{J}\left\langle\hat{ \mathsf{n}}\right\rangle-1+\frac{\left\langle\hat{\delta}_{n}\hat{\delta}_{n} \right\rangle}{\left\langle\hat{\mathsf{n}}\right\rangle} \tag{24}\]  
occur when \(\widetilde{\left\langle 1\right\rangle}=R\). The full set of equations of motion can be found in Appendix B, with an analytical solution for the \(F=0\) case. Expectation values such as Equation (22) are given in Appendix C. From a different starting point, coupled equations for the evolution of density, phase and their variances have also been used in [67].