# Gaussian Quantum Trajectories for the Variational Simulation of Open Quantum-Optical Systems
## Appendix B Full Equations for \(N\Theta\)-Gaussian Trajectories  
For the photon-counting unraveling, expectation values evolve as:  
\[\partial_{t}\left\langle\hat{\mathsf{n}}\right\rangle = -\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle+2\operatorname{Im}\left[FC_{1}\right] \tag{12}\] \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\] \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\] \[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\] \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{ \theta}\right\rangle_{\text{sym}} = -U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n}\right\rangle +2\operatorname{Im}\left[FC_{6}\right]-\operatorname{Re}\left[FC_{7}\right]\] \[\partial_{t}\widehat{\left\langle\hat{1}\right\rangle} = -\gamma\left\langle\hat{\mathsf{n}}\right\rangle\widehat{\left\langle \hat{1}\right\rangle},\]  
where correlators \(C\) are given with an expansion in \(\left\langle\hat{\mathsf{n}}\right\rangle\) in Appendix C.  
Jumps when \(\widehat{\left\langle 1\right\rangle}=R\) are given by:  
\[\left\langle\hat{\mathsf{n}}\right\rangle \xrightarrow{J}\left\langle\hat{\mathsf{n}}\right\rangle-1+ \frac{\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n}\right\rangle}{ \left\langle\hat{\mathsf{n}}\right\rangle} \tag{13}\] \[\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n}\right\rangle \xrightarrow{J}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle\left(1-\frac{\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n }\right\rangle}{\left\langle\hat{\mathsf{n}}\right\rangle^{2}}\right)\] \[\left\langle\hat{\mathsf{s}}\right\rangle \xrightarrow{J}\left\langle\hat{\mathsf{s}}\right\rangle+\frac{ \left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle_{\text{ sym}}}{\left\langle\hat{\mathsf{n}}\right\rangle^{2}}\] \[\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}} \xrightarrow{J}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{ \theta}\right\rangle_{\text{sym}}\left(1-\frac{\left\langle\hat{\mathsf{s}}_{n }\hat{\mathsf{s}}_{n}\right\rangle}{\left\langle\hat{\mathsf{n}}\right\rangle^{2 }}\right).\]In the \(F=0\) case, the evolution (A4) can be solved exactly to:  
\[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]0 (A6) \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]1 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]2 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]3 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]4 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]5  
By equating \(\widetilde{\left\langle 1\right\rangle}(t_{j})=R\), we find the time to the next jump to be:  
\[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]6 (A7)  
In the heterodyne unraveling, equations for the evolution are given by:  
\[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]7 (A8) \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]8 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{n} \right\rangle = 4\operatorname{Im}\left[FC_{2}\right]-2\operatorname{Im}\left[FC _{1}\right]\]9 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]0 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]1 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]2 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]3 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]4 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]5 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]6 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]7 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]8 \[\partial_{t}\left\langle\hat{\mathsf{s}}\right\rangle = \left(\Delta+\frac{U}{2}\right)-U\left\langle\hat{\mathsf{n}} \right\rangle-\gamma\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta} \right\rangle_{\text{sym}}-\operatorname{Re}\left[FC_{3}\right]\]9  
\[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\]0 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\]1 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\]2 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\]3 \[\partial_{t}\left\langle\hat{\mathsf{s}}_{\theta}\hat{\mathsf{ s}}_{\theta}\right\rangle_{\text{sym}} = -2U\left\langle\hat{\mathsf{s}}_{n}\hat{\mathsf{s}}_{\theta}\right\rangle _{\text{sym}}-2\operatorname{Re}\left[FC_{4}\right]+\frac{1}{2}\operatorname{ Im}\left[FC_{5}\right]\]4  
where for the expansion up to order \(\left\langle\hat{n}\right\rangle^{-1}\) between the brackets, all Wick contractions are taken into account (re-summing over all orders in \(\hat{\delta}_{\theta}\)), and for order \(\left\langle\hat{n}\right\rangle^{-2}\), higher order correlators between the prefactor and exponential are neglected, similar to [84].