Gaussian Quantum Trajectories for the Variational Simulation of Open Quantum-Optical Systems

We construct a class of variational methods for the study of open quantum systems based on Gaussian ansatzes for the quantum trajectory formalism. Gaussianity in the conjugate position and momentum quadratures is distinguished from Gaussianity in density and phase. We apply these methods to a driven-dissipative Kerr cavity where we study dephasing and the stationary states throughout the bistability regime. Computational cost proves to be similar to the truncated Wigner (TWA) method, with at most quadratic scaling in system size. Meanwhile, strong correspondence with the numerically exact trajectory description is maintained so that these methods contain more information on the ensemble constitution than TWA and can be more robust.


INTRODUCTION
Many-body systems of interacting photons have come under intense investigation over the last years, with both circuit QED and semiconductor heterostructure based systems [1][2][3][4][5][6]. The main difference with traditional many-body systems, both in hard and synthetic condensed matter, is the fact that the photon lifetime is typically shorter than the time at which the system dynamics develop. To compensate for the losses, photonic systems have to be continuously driven, yielding a steady state that is a balance between driving and dissipation. Even though the basic theoretical framework for the description of open systems is well understood [7,8], the exponential complexity of the many-body problem requires the development of practical approximate techniques to tackle systems that consist of many excitations of multiple modes, stirring the need for methods that unite the approaches of these different fields [9].
In order to study the time-evolution of an open quantum system, there are two distinct approaches [8]. The first one is the study of the ensemble as a whole. For a Markovian system this translates into the Lindblad master equation for the density matrixρ, whereρ contains D 2 elements for Hilbert-space dimension D. The second approach is the trajectory approach: here only a pure wave-function |ψ of D elements is evolved at a time. Through stochastic decoherence, the environment effectively introduces additional classical noise that complements a non-hermitian Schrödinger evolution of the system, and the ensemble evolution is obtained by averaging over many trajectory evolutions. In addition to the reduction of complexity from D 2 to D, another advantage of the trajectory approach is that the full statistics of detector clicks are obtained.
After preliminary work by Davies [10], this quantum trajectory formalism was developed towards its current shape [7,11] by a number of different groups [12][13][14][15]. It is also known under a variety of other names: Monte- * wouter.verstraelen@uantwerpen.be Carlo wave function method, quantum jump method or Stochastic Schrödinger equation. The act of averaging over many such trajectories to obtain an ensemble description is often called the Stochastic simulation Method. The formalism behind quantum-trajectories is tightly bound to quantum metrology [16]. In this context, is important to note that quantum trajectories are not a single method, but rather a class of methods, depending on the so-called unraveling which correspond to some (hypothetical) measurement protocol taking place on the environment. For quantum optical systems, the most common protocols are photon-counting, resulting in piecewise-deterministic processes on the one hand and homodyne/(far-detuned-)heterodyne detection resulting in Wiener processes on the other hand [8]. The heterodyne case has proven to be equivalent to the stochastic collapse model of quantum state diffusion [8].
Although analytical effort can sometimes provide a serious reduction [17,18], computational complexity typically remains a limiting factor in exact numerical study of open systems. The aforementioned quantum trajectory method, which allows for a description in the Hilbert space of pure states H instead of the superspace of density matrices H 2 , can reduce the memory-cost for individual simulations to the level of closed-system problems, at the expense of the need to average over many individual Monte Carlo realizations. Still, as in the closed system case, many-body systems are usually too complex for exact computation and one needs to use approximate methods instead. These can either occur at the level of the Master equation or at the level of individual trajectories [9]. We are interested in a method that is variational and provides a description on the level of trajectories. Some current variational (TDVP, Time-Dependent Variational Principle [19]) approaches to the description of open systems aim to describe the system at the level of the Master equation [20][21][22][23][24], with ansatzes including including Gutzwiller Density Matrix- [25][26][27], Matrix-Product State [28,29] and Matrix-Product Operator [30][31][32][33] methods. Variational descriptions on the trajectory level have so far focused towards 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 [34][35][36] or based on t-DMRG [37][38][39].
Somewhat surprisingly, except for some works on the back-reaction of measurements on quantum many-body systems [40,41] , the class of Gaussian states [42] 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 [43,44], and with the important application of time-dependent BCS/BdG [45,46] theory. It is the purpose of this work to develop the formalism of a Gaussian variational description of quantum trajectories.
This description is expected to work well in the regime of many particles and weak interactions. In this regime, a different stochastic technique, the so-called truncated Wigner approximation (TWA) is widely used [47][48][49][50][51], which fits in a broader context of phase-space methods [42,47]. Here, the starting point is the Fokker-Planck equation for the Wigner function of the system. If one neglects derivatives of order higher than two, this Wigner function can be sampled through Feynman-Kac onto stochastic differential equations (an alternative formulation for open systems is given in [52]). As opposed to the aforementioned trajectory approach, single samples of TWA do not correspond to a physical state but rather sample the Wigner phase-space. Though TWA can yield very good results, it is not always well-controlled and may predict unphysical behavior [53].
With this work, we wish to close the gap between the robustness of exact quantum trajectory methods and the numerical efficiency of the TWA method. It is organized as follows: in the next section II, we refer to the quantum trajectory method, how it affects expectation values and the meaning of a Gaussian ansatz. In III we translate this idea to the states that are commonly known as '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 [54]. 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 IV 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 are N Θ-Gaussian methods, meaning that the state is Gaussian in density and phase. In section V we will adress computational issues and compare performance with TWA and exact trajectory methods. Finally, in section VI 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.

A. Quantum trajectories for expectation values
For the remainder of this work, we assume for simplicity systems where the interaction with the environment is entirely characterized by their Markovian dissipation (photon-leaking towards empty space) at rate γ. At the level of the master equation for the ensemble density matrixρ, such a system with HamiltonianĤ is described by withâ(â † ) the photonic annihilation(creation) operators.

Photon-counting unraveling
According to the quantum trajectory method, this master equation (1) is unraveled into the evolution of stochastic wavefunctions by performing continuous measurement on the environment. The simplest unraveling is photon-counting (PC ) detection: in between detectorclicks, the wavefunction is propagated as where the tilde denotes that the norm of the wavefunction is not conserved. The detection of a photon corresponds to a discrete jump when the norm has decreased to |ψ 2 = R, where the random number R has a uniform distribution on the interval [0, 1].
In terms of expectation values, equation (2) corresponds to the evolution By choosingÔ = 1, we see that the evolution for the normalized expectation value Ô = Ô / 1 is given by Meanwhile, according to (3), a jump affects the expectation value Ô by equations of the form (5) and (6) together then in principle provide a complete description of the evolution of expectation values for a trajectory under photon-counting unraveling.

Homo-and heterodyne unravelings
Apart from photon-counting, (generalized) homodyne detection can be performed, where the emitted light is interfered with a classical reference signal (local oscillator) β = |β|e iω LO t , resulting in the measurement of the quadrature variables X and P . Homodyne measurement corresponds to jump operatorsĴ = √ γ(â +β).
One readily finds that by inverting this relation towardŝ a =Ĵ √ γ − β and substituting in (1), again an equation in the Lindblad form is retrieved when absorbing a contribution √ γ Im[β * Ĵ ] in the Hamiltonian. Therefore, homodyne detection is equally valid as an unraveling of the Master equation. In practice, because β is macroscopic, a diffusion approximation of the jumps is justified [8], yielding the Itô equation for the unnormalized wavefunction |ψ . Here we have divided the decay rate in two channels γ = γ X + γ P for which either X or P are monitored. That is, if β is inphase with the cavity field, there is only one independent measurement record and we speak of 'true' homodyne detection, for example if only X is measured (Hom. (X)), γ X = γ and γ P = 0. If β is far-detuned from the cavity field on the other hand, both quadratures are equally measured such that γ X = γ P = γ 2 . This latter detection scheme is known as heterodyne detection (Het.). From the next section on, we will assume this heterodyne case unless indicated otherwise. The complex Wiener noise in Eq. (7) has a Gaussian distribution with zero mean and variance dZ = γ X γ dW X + i γ P γ dW P with dW 2 X = dW 2 P = dt and dW X dW P = 0, meaning that |dZ| 2 = dt and, in the heterodyne case, dZ 2 = dZ * 2 = 0.
Analogously to the deterministic part of photoncounting, (7) can be translated to the evolution of nor-malized expectation values, yielding where we have introduced the generic fluctuation

B. Closing at Gaussian level
Equations (5) or (8) will generally introduce an infinite hierarchy of correlation functions. This is a similar situation to classical mechanics, where the Liouville equation introduces the BBGKY hierarchy, that must be truncated at some point. At the level of evolution for the full ensemble (the master equation), a systematic derivation of proper truncation schemes (also known as cumulant expansions) for driven-dissipative quantum systems has been performed in [55]. Generally, one expects the resulting predictions to converge as more correlation functions of higher order are taken into account. Nevertheless, it is only at the mean-field or at the Gaussian level that it is possible to close the hierarchy instead of truncating, ensuring that the state remains physical by construction. When closing (5), (8) at mean-field level, the dynamics are trivial and coincide with the mean-field solution of the master equation. At the Gaussian level on the other hand, the closing of the equations can be performed by applying Wick's theorem and does reflect the stochasticity of the unravelings.
As an ansatz, we thus approximate the state to be contained in a Gaussian subspace of the Hilbert space. These Gaussian states are entirely characterized by their first and second moments, keeping the amount of independent variables of which the evolution needs to be studied limited. This provides a large reduction of complexity from the whole Hilbert-space, which scales exponentially with system size (the number of modes). Because of the presence of the second moments, a Gaussian ansatz is by nature one order higher than the mean-field approach that corresponds to a coherent ansatz. Formally then, the Gaussian methods are of the same order as Bogoliubov theory. Nevertheless, we can expect these Gaussian trajectory methods to provide a more refined description of the whole state because the second-order expansion is done at the level of individual trajectories instead of on the level of the ensemble. In the next section, we will illustrate this with the example of bistability in a Kerr nonlinear cavity. Where the Bogoliubov theory is only able to describe the fluctuations around one of the two stable branches, the variational Gaussian method is also capable of describing the switching between the branches.

III. KERR-BISTABILITY AND THE XP -GAUSSIAN METHODS
As a first example, we look at a driven-dissipative cavity with Kerr non-linearity, described by the Hamiltonian where ∆ the cavity-laser detuning, U the photon-photon interaction and F the classical laser amplitude. Furthermore, photons leak to the vacuum with rate γ as described in section II. Exact solutions for the expectation values â † mâ n of the stationary state, as well as a semi-classical model, were calculated in [54]. Interestingly, at the classical level, this system features a bistability regime where, given all other parameters, two stable solutions for the density exist.
Such a Kerr-cavity provides a model for polariton con-densates [47], for which it has been shown that trajectory methods can provide an adequate description [56,57]. Also coupled arrays of these cavities are an object of current interest [58] and it is in these systems that the TWA method has been proven insufficient [53].
As it is the most straightforward ansatz corresponding with the Gaussian approximation of section II, we assume the the state to be Gaussian in the quadrature operatorsX andP or, equivalently,â andâ † , i.e. the states that are colloquially known as Gaussian states [59]. Such an XP -Gaussian state is characterized entirely by its mean-field value α = â and its two-point correlation functions ââ and â †â or, equivalently, α, δδ and δ †δ whereδ =δâ. The connected correlation functions equal δδ = ââ − α 2 and δ †δ = â †â − |α| 2 .
From (5) and after applying Wick's theorem to close the set of equations, we obtain for photon-counting the deterministic evolution This deterministic evolution is propagated until the norm, evolving through (4), becomes 1 = R. Then, according to (6), jumps occur as In the heterodyne unraveling on the other hand, the evolution of the mean-field (8) is given by Generally, one would expect similar stochastic equations for δδ and δ †δ . Remarkably, they are entirely deterministic and, moreover, are identical to the deterministic part of photon-counting equations (14), (15). It should be noted, however, that the the equations of motion are different for more general homodyne detection schemes.
For all unravelings, as a slightly more efficient and stable alternative to evolving δ †δ along explicitly trough (15), one can obtain it from δδ by asserting that the state remains pure, corresponding to the condition [59] δ †δ + δ †δ 2 = δδ 2 .
For an initially pure state, relation (17) remains exactly satisfied through heterodyne measurement as well as through the deterministic evolution between photon jumps, although the purity briefly decreases at a jump after which it tends to relax back to 1. Relation (17) remains fulfilled at a jump up to order |α| −5 so that only in systems where there is a significant jumprate at zero density it is better to evolve δ †δ explicitly with Eq. (15). In Fig. 1, we show a comparison between single trajectories that were obtained with the numerically exact evolution of the wave function and the Gaussian variational ansatz where for the jumps in both simulations identical random numbers were used. Initially there is a very strong correspondence with the XP -Gaussian variational method while at later evolution times an offset in time emerges and the good correspondence is lost. Qualitatively though, the behavior remains similar so that approximate correspondence in the weak sense (for the whole ensemble) can remain fulfilled.
To verify this, in figure 2 densities and density correlations g (2) = â †â †ââ / â †â 2 are plotted as function of FIG. 2. Stationary ensemble expectations of photon number (left) and second-order correlation (right) as function of F through the bistability regime for parameters U/γ = 0.05, ∆/γ = 1. Results were obtained by averaging over 10 4 samples after t = 100γ −1 of evolution from the vacuum.
F throughout the bistability regime of the Kerr-model. We see that the XP -Gaussian methods can provide reasonable to very good predictions of the exact correlation functions from [54]. In this regime of relatively low density these XP -Gaussian methods are mostly outperformed by the TWA method, though. It is also seen that that the predictions on the statistics of the ensemble as a whole depend (weakly) on the choice of unraveling, whereas this is independent for exact trajectories. We should note that, because of the low density, simulation with exact trajectories can also still easily be performed. It is the opposite, high-density, limit where Gaussian ansatzes will provide a better description as well as where exact trajectories become computationally unfeasible.

IV. PHASE DIFFUSION AND THE
N Θ-GAUSSIAN METHOD A. Phase space evolution As another example, we look at the time-evolution of a freely evolving Kerr cavity: we envision an initial state present in system (9) where at t = 0 the pump is turned off, i.e. F = F on Θ(−t). Without the pump, U(1) symmetry is restored, allowing the phase to diffuse freely.
The challenging nature of this problem for our Gaussian variational ansatz can be appreciated from the Wigner distributions of a single realization after some time, shown in Fig. 3. Panel (a) shows that the phase space distribution has almost spread out over a full circle, implying the loss of phase coherence. The difference between the left and right hand panels, obtained with photon counting and heterodyne detection respectively, shows how the phase space distribution is kept more concentrated under heterodyne measurement as compared to the photon counting. This is expected, because only the heterodyne measurement gives phase information.
In this example, we clearly see the importance of the unraveling on the applicability of a Gaussian approximation. Unfortunately for the parameters of Fig. 3, even for It is clear to see that these states are not XP -Gaussian, firstly because they are bended and secondly because the W -function exhibits Fock-like negative parts. Note that the intra-sample variance Var1 of the phase is smaller for a heterodyne sample, its inter-sample Var2 variance is larger. Therefore, heterodyne detection is slightly less problematic for the XP −Gaussian states.
the heterodyne measurement, a Gaussian approximation is very crude. The applicability of the Gaussian method is directly related to the intra sample variance of the phase. In general, in a quantum trajectory simulation, the variance of an observable can be written as where the intra and inter trajectory variances are respectively defined as [8] Var where index α labels the trajectories. Only Var(Ô) can be measured without unraveling the dynamics and it is the only one that is accessible within the master equation description (and hence the TWA).
The parameter regime where the width of the phase space distribution is expected to become small can be found by requiring that the phase diffusion rate is much smaller than the rate at which phase information is obtained: U â †â γ â †â . For the parameters of Fig.  3, we have U/( â †â γ) = 0. 1 1, but still the phase space distribution cannot be accurately approximated by a Gaussian. Mean (left panel) and inner-sample-, betweensample-and total variances (right panel) of the X-quadrature obtained by the exact trajectory methods and TWA as well as XP -Gaussian and N Θ-Gaussian variational methods as function of time after free evolution from a coherent initial state |α = |10 (same parameters as Fig. 3). Averages are taken over 1e3 samples (1e4 for TWA). It is clear that the N Θ variational method provides an accurate description of the true dynamics, both on the level of the whole ensemble as on the level of its constitution in pure states, whereas the latter cannot be obtained from the TWA method. Very similar behavior is present regarding P and its variances as well as for the X, P −covariances.

B. N Θ-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 [60]. Nevertheless, we will work with the Dirac definitionâ =: e iθ √n , whereθ defines the phase andn is the familiar particle number operator. The practical use of this phase is paramount in quantum hydrodynamics [61]. 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: n,θ = i.
Assuming the state to be Gaussian inn andθ, the independent (real) expectation values to take into account are n , δ nδn = nn − n 2 , θ , δ θδθ = θθ − θ 2 and δ nδθ sym = nθ /2 + θn /2 − n θ . From these correlation functions, expectation values of products of annihilation and creation operators can be approximately computed by expanding the exponential and using Wick's theorem. For example, 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 II. For example, for the deterministic evolution under photon-counting we obtain for the density which is the equivalent of (17). For a state that is initially pure, relation (26) again remains satisfied up to order n −5 after a jump. Equations of motion for heterodyne measurement of density-phase Gaussian states are also given in appendix A, for example for the density, d n = 2 Im F √n e -iθ − γ n dt is obtained. On figure 4 (a) the expected evolution of quadrature variableX is shown. The right panels show the intrasample variance Var 1 , inter-sample variance Var 2 and total variances for the evolution of a state that was originally coherent. We see here that the XP -Gaussian methods for both the photon-counting and the heterodyne unraveling strongly underestimate the phase diffusion, the description of which is worse for the photon-counting description.
Only in regimes where the importance of losses is much higher than the importance of dephasing can XP -Gaussian states provide a suitable description. On the other hand, we see that the N Θ-methods are able to capture the phase diffusion on the level of the ensemble as well as TWA. What distinguishes the former from the latter however, is that the N Θ-Gaussian method is able to show the composition of the ensemble: it maintains information of individual trajectories, which is lost in TWA.

V. COMPUTATIONAL ASPECTS
Because of the few simple update rules (five real variables-four if purity relations are used), the computational cost for evolving a single variational trajectory is of the same order as the cost for a single TWA sample (two real variables) given the same set of parameters. For an exact trajectory D = 2N levels real variables must be evolved for the same single mode, where N levels is the amount of occupation levels considered. If the dimension of configuration space d increases, the dimension of Hilbert space correspondingly grows as D ∝ e d , so that large simulations easily become limited by computational load. The dimension of the corresponding phase space on the other hand only grows ∝ d, making it far more efficient for large systems. For our variational methods it is the amount of distinct Gaussian correlation functions that is the relevant quantity of complexity and this scales in principle as ∝ d 2 . In many large systems, only correlations between neighboring sites are important and other ones can be neglected. If this is the case, scaling ∝ d is retrieved, as in Gutzwiller-or tensor-network ansatzes.
So far, we have only discussed the evolution of a single sample. As Gaussian trajectories, like exact trajectories, have a finite spread in Wigner phase-space, they contain more information than a TWA-sample which is just a point. Therefore, it can be expected that less samples are required to obtain an accurate description of the ensemble dynamics. We can estimate the difference in statistics as follows: take a trajectory method (exact or variational) of which N traj independent trajectories are evolved. For the trajectory methods, the statistical uncertainty [62] is whereas the analogous uncertainty for TWA is The ratio of N traj and N TWA needed for the same precision can be found by equating (28) and (29) to be In this sense, variational trajectories may even outperform TWA computationally. As an example, on figure 5 the same process as on figure 4 is simulated with only ten samples of the exact/variational trajectories and ten samples of TWA. It is clear that TWA provides a bad estimate with such a low number of samples. The performance of trajectory methods with regard to the minimal  figure 4 (corresponding legend) simulated with ten samples of exact/variational trajectories as well as TWA. It is clear that trajectory methods are more tolerant for low statistics than TWA. Corresponding to criterion (30), we see that photon-counting is still more tolerant than heterodyne detection, corresponding with the fact that Var1(X) Var2(X) for this photon-counting unraveling, as opposed to heterodyne. Note that the initial state already has a deviation from the true value for TWA but not for the trajectory methods: criterion (30) for the needed ratio of samples diverges there.
amount of samples, though observable-and unravelingdependent, is better. The amount of improvement depends onÔ and the unraveling: as seen from Eq. (30), the photon counting unraveling performs best because Var 1 (X) Var 2 (X) as can be appreciated from figure  4.
When it is a phase-space distribution itself that one is interested in, it can again be challenging to compute it from an exact density matrix or wavefunction, when the particle number is not low [63]. If a state is Gaussian, particularly in the XP -sense, distributions for individual samples are given by straightforward Gaussian formulas [42] that can be added up for the distribution of the whole ensemble. The TWA method is by construction surely also very well suited to plot the Wigner-quasi-probability function, but binning as a histogram is necesarry there, meaning that again a much larger amount of samples is required for a result of high resolution.
As a final note, all simulations discussed have been performed by trajectories representing pure states. It is also possible to work with trajectories using mixed states, either corresponding to initial classical uncertainty or to imperfect measurements [8]. The exact solution of these would require a density matrix, undoing the computational advantage of the stochastic simulation method. For the variational trajectories, the only difference between evolving pure or mixed states is whether constraints (17) and (27) are valid or not. Evolving mixed states increases Var 1 (Ô) with respect to Var 2 (Ô) so that less statistics are required according to (30), at the expense of resolution of the individual trajectories . This would come down to a crossover from pure trajectories to a description on the level of the master equation. Whether such a description is accurate will depend on the system dynamics. For the example of the bistabil-ity, if there is not sufficient information on the branch that the system is in, the Gaussian approximation will fail to accurately describe the state, that evolves toward a mixture of the system in the lower and upper branch.

VI. CONCLUSIONS AND OUTLOOK
We have shown how Gaussian variational quantum trajectory methods can provide a computationally efficient description of open quantum systems. We have explicitly derived the dynamics for both XP -and N Θ-Gaussian states, and applied them as an example to a drivendissipative cavity. XP -Gaussian states are always welldefined, though they may be too rigid to describe states with much phase-diffusion. N Θ-states on the other hand, exist only by approximation, but often provide a very good description for the true state as long as the density is sufficiently high. Computationally, the cost of a Gaussian trajectory scales similar to TWA, but typically less samples are needed. XP -Gaussian methods, like TWA, are limited to the semiclassical regime where the Wigner function is always positive. N Θ-Gaussian methods on the other hand can describe the interference patterns similar to Fock-states (XP -Gaussian states are the only pure states with an entirely positive Wigner function [59]). As these variational methods behave more controlled than TWA, we expect them to keep predicting accurate results to a broader class of problems such as coupled cavities. As we have shown, accuracy can be strongly improved by proper choice of the ansatz. Also similar ansatzes, that can for example be Gaussian in other variables [64], may be suitable dependent on the problem. As for exact trajectories, the choice of unraveling greatly determines the amount of samples needed for a proper description. In addition, we have seen that this choice of unraveling can also influence the accuracy of the variational method. Besides the extension to larger systems, the application to systems with multiple Markovian jump processes is straightforward. Extensions towards more general non-Markovian noise are less clear-cut, but may be feasible through e.g. a doubled Hilbert space framework [8].

ACKNOWLEDGMENTS
We thank Wim Casteels for his important initial contributions to this project; and acknowledge discussion with M. Van Regemortel and D. Huybrechts. This work was financially supported by the FWO Odysseus program. Part of the computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation -Flanders (FWO) and the Flemish Government department EWI.
In the heterodyne unraveling, equations for the evolution are given by