IInterference of Quantum Trajectories
Brecht Donvil and Paolo Muratore-Ginanneschi
University of Helsinki, Department of Mathematics and Statistics P.O. Box 68 FIN-00014, Helsinki, Finland ∗ We show that general time-local quantum master equations admit an unravelling in quantum trajectories withjumps. The sufficient condition is to weigh state vector Monte Carlo averages by a probability pseudo-measurewhich we call the ”influence martingale”. The influence martingale satisfies a d stochastic differential equationenslaved to the ones governing the quantum trajectories. Our interpretation is that the influence martingale mod-els interference effects between distinct realizations of the quantum trajectories at strong system-environmentcoupling. If the master equation generates a completely positive dynamical map, there is no interference. Insuch a case the influence martingale becomes positive definite and the Monte Carlo average straightforwardlyreduces to the well known unravelling of completely positive divisible dynamical maps. PACS numbers: 03.65.Yz, 42.50.Lc
Background and motivation
Actual quantum systems areopen: they unavoidably interact, even if slightly, with theirsurrounding environment [1]. A useful phenomenological ap-proach is to conceptualize the interaction as a generalizedmeasurement performed by the environment onto the system[2, 3]. As a consequence, the state vector of an open systemfollows stochastic trajectories in the Hilbert space of its styl-ized isolated counterpart. These trajectories are characterizedby sudden transitions, quantum jumps. Since the experimen-tal breakthroughs [4, 5] quantum jumps have been observed inmany atomic and solid-state single quantum systems see e.g.[1] for an overview. Recent experimental developments evensupport the possibility to identify precursors of the imminentoccurrence of a jump [6].From the theoretical slant, the connection between quantumtrajectories and the master equation formalism is very well un-derstood under general but non all comprehensive conditions.In particular (but not exclusively, see end of this introduction),a rigorous mathematical theory exists if the state operator ofthe open system obeys a completely positive divisible dynam-ics [7–11]. Remarkably, similar results have also been inde-pendently derived in the study of the foundations of quantummeasurement theory [12–15].The state operator of an open system evolves according toa completely positive linear dynamical map ρ t = Φ t ( ρ ) ifand only if Φ t is specified by the partial trace over the envi-ronment of the unitary dynamical map U t ρ ⊗ σ U † t [16].Here U t is the fundamental solution of the Schr¨odinger equa-tion of the closed system-environment evolution from time to t . Thus, the initial state at time must be a tensor productof state operators acting respectively in the Hilbert space ofthe system ( ρ ) and the environment ( σ ) .If a dynamical map Φ t is invertible at any time, then it isalso divisible [17] Φ t s = Φ t Φ − s (1)In such a case and, assuming a smooth time dependence, Φ t s satisfies a time local master equation Φ t + ε s − Φ t s = ( Φ t + ε t − Φ t s = G t ( Φ t s ) ε + o ( ε ) with generator G t . In the weak coupling van Hove scalinglimit [2, 18, 19] the generator G t takes the form G t ( ρ t ) = − ı [H , ρ t ] + D( ρ t ) (2a) D( ρ t ) = (cid:76) (cid:88) (cid:96) =1 Γ (cid:96),t (cid:16)(cid:104) L (cid:96) , ρ t L † (cid:96) (cid:105) + (cid:104) L (cid:96) ρ t , L † (cid:96) (cid:105)(cid:17) . (2b)The limit procedure determines the explicit expression of theeffective Hamiltonian H (eventually time dependent) and ofthe Lindblad ”dissipative” operators { L (cid:96) } (cid:76) (cid:96) =1 ( (cid:76) arbitrary)based on the unitary dynamics of the microscopic system-environment model. The Γ (cid:96),t ’s are time dependent functionsweighing the Lindblad operators. More generally, Lindblad[20] and Gorini, Kossakowski and Sudarshan [21] (see also[22] for an insightful brief overview) proved that a masterequation of the form (2) is the generator of a completely pos-itive map if and only if for all times t Γ (cid:96),t ≥ (cid:96) = 1 , . . . , (cid:76) . (3)The conditions (3) specify a completely positive divisible evo-lution. The connection with quantum trajectory theory is thatif (3) holds true, the solution of (2) can be computed as aMonte Carlo average, denoted here by E , over a stochasticstate vector evolution in the system Hilbert space ρ t = E ψ t ψ † t (4)At variance with classical master equations, (2) does not fixthe statistics of the stochastic process { ψ t } t ≥ . If quantumjumps are modelled by Poisson processes, ψ t obeys an Itˆostochastic Schr¨odinger equation of the form [3, 10, 23, 24] d ψ t = d t f t + (cid:76) (cid:88) (cid:96) =1 d ν (cid:96),t (cid:18) L (cid:96) ψ t (cid:107) L (cid:96) ψ t (cid:107) − ψ t (cid:19) (5a) ψ = z (5b)with z † z = 1 . The { ν (cid:96),t } (cid:76) (cid:96) =1 are a collection of independentPoisson processes (see e.g. [3]) and their increments satisfy d ν (cid:96),t d ν (cid:107) ,t = δ (cid:96), (cid:107) d ν (cid:96),t (cid:96), (cid:107) = 1 , . . . , (cid:76) (6a) E (cid:0) d ν (cid:96),t (cid:12)(cid:12) ψ t , ¯ ψ t (cid:1) = (cid:114) (cid:96),t (cid:107) L (cid:96) ψ t (cid:107) d t (cid:96) = 1 , . . . , (cid:76) (6b) a r X i v : . [ qu a n t - ph ] F e b Condition (3) allows us to identify the rates of the Poissonprocess with the weights of the Lindblad operators in (2) (cid:114) (cid:96),t = Γ (cid:96),t (cid:96) = 1 , . . . , (cid:76) (7)Thus, the choice f t = − ı H ψ t − (cid:76) (cid:88) (cid:96) =1 Γ (cid:96),t L † (cid:96) L (cid:96) − (cid:107) L (cid:96) ψ t (cid:107) ψ t (8)and a straightforward application of Itˆo calculus for Poissonnoise [25] to compute the differential of (4) recover the masterequation (2).These beautiful results, however, do not encompass all timelocal master equations. The inverse of a completely posi-tive map is completely positive only if the evolution is uni-tary [19]. Therefore, a completely positive map is not ingeneral completely positive divisible. And if (3) does nothold true, the bridge relation (7) between quantum trajecto-ries and master equation cannot be satisfied. Explicit exam-ples are time local master equations obtained from Gaussiansystem-environment models of Bosons e.g. [26] or Fermionse.g. [27, 28], other exactly integrable models of system-environment interactions e.g. [29], and more generally masterequations obtained by means of time convolutionless pertur-bation theory [2]. In other words, the exact or approximatetime local evolution law of Φ t s for generic s may be oftenmore accessible than identifying the time s = 0 when systemand environment are exactly represented by the tensor productof two independent state operators.In order to introduce quantum trajectories when (3) does nothold, [30, 31] interpreted negative values of the Γ (cid:96),t as ”cor-relations” between the components of the state vector of thereduced system and those of an auxiliary system in a doubledHilbert space. More recently, [32, 33] (see also [34]) derived aquantum trajectory description at the level of the system stateoperator . Such description has the advantage to avoid compu-tationally taxing extensions of the Hilbert space. The resultsof [32, 33] require, however, to keep track of memory effectsbetween distinct realizations of the stochastic state operator.Here we show that a re-definition of (4) which accountsfor interference, constructive or destructive, between quantumtrajectories permits to unravel (2) when (3) does not hold.Before turning to the description of our result, we want toemphasize that the notion of quantum trajectory is not nec-essarily related to time local master equations. Diosi, Strunzand Gisin [35, 36] derived the unravelling of non-necessarilyinvertible dynamical maps occasioned by the tracing out anenvironment consisting of a bath of boson oscillators. Theirresults have been further refined [37, 38], including fermionbaths [39] and also exploring its significance in connection tothe modal interpretation of quantum mechanics [40, 41]. Fi-nally, a quantum trajectory can be associated to the continuousmonitoring of an output field commuting with itself at differ-ent times [42]. The reduced system state operator can thenbe unravelled in terms of a complete set of states of the reser-voir corresponding to a continuous measurement of the outputfield. Main result.
The idea of the influence martingale methodis to compute the state operator as a weighed Monte Carloaverage over the Poisson measure: ρ t = E (cid:16) µ t ψ t ψ † t (cid:17) (9)The weight { µ t } t ≥ is a martingale process satisfying thescalar Itˆo stochastic differential equation d µ t = µ t (cid:76) (cid:88) (cid:96) =1 (cid:103) (cid:96),t (cid:16) d ν (cid:96),t − (cid:114) (cid:96),t (cid:107) L (cid:96) ψ t (cid:107) d t (cid:17) (10a) µ = 1 (10b)driven by the same Poisson processes (6) of (5) and by a col-lection { (cid:103) (cid:96),t } (cid:76) (cid:96) =1 of pure functions of time whose explicitform we determine below. The process { µ t } t ≥ is a (local)martingale because the Poisson processes are compensated bytheir conditional expectation so that at all times E(d µ t | ψ t , ¯ ψ t ) = 0 . More precisely, (10a) by itself only describes a local martin-gale (see e.g. [25]). A local martingale becomes a strict mar-tingale if the integrability condition
E sup t | µ t | < ∞ holds.On physics grounds, we expect the integrability condition tohold whenever a time local master equation of the form (2)is well defined. We thus identify in what follows the no-tions of local and strict martingale. In that case, (10b) implies E µ t = 1 for all t . Under these assumptions, a straightforwardapplication of stochastic calculus yields dd t E( µ t ψ t ψ † t ) = − ı [H , ρ t ] − (cid:76) (cid:88) (cid:96) =1 Γ (cid:96),t L † (cid:96) L (cid:96) ρ t + ρ t L † (cid:96) L (cid:96) (cid:76) (cid:88) (cid:96) =1 (cid:16) ˜Γ (cid:96),t L (cid:96) ρ t L † (cid:96) − (˜Γ (cid:96),t − Γ (cid:96),t ) E( (cid:107) L (cid:96) ψ t (cid:107) µ t ψ t ψ † t ) (cid:17) with ˜Γ (cid:96),t = (1 + (cid:103) (cid:96),t ) (cid:114) (cid:96),t . Upon contrasting this result with(2) we obtain a new bridge relation between master equationand quantum trajectories Γ (cid:96),t = ˜Γ (cid:96),t ≡ (1 + (cid:103) (cid:96),t ) (cid:114) (cid:96),t (11)Until now, we did not impose any conditions on the (cid:103) (cid:96),t ’s.Thus, drawing from [32] we may satisfy (11) by requiring (cid:114) (cid:96),t = | Γ (cid:96),t | (12a) (cid:103) (cid:96),t = sign(Γ (cid:96),t ) − (12b)With these choices, the influence martingale method imme-diately reduces to a standard unravelling ( µ t = 1 for all t )when (3) holds true. Furthermore (12) exhibit that µ t dependsupon ψ t without exerting any feedback. In this sense (10) is enslaved to the stochastic Schr¨odinger equation (5).The question then arises about the physical interpretationof the Hilbert space vector valued process { ψ t } t ≥ whenthe rates of the Poisson processes (6) do not satisfy (3). Thesquared norm satisfies the stochastic differential equation d (cid:16) (cid:107) ψ t (cid:107) (cid:17) = (cid:76) (cid:88) (cid:96) =1 (cid:0) d ν (cid:96),t − Γ (cid:96),t (cid:107) L (cid:96) ψ t (cid:107) d t (cid:1) (cid:16) − (cid:107) ψ t (cid:107) (cid:17) which is not a martingale unless (3) holds true. Nevertheless,the (hyper)-Bloch sphere (i.e. the manifold (cid:107) ψ t (cid:107) = 1 ) ispreserved by the dynamics. Thus, we can interpret ψ t as statevector at time t and then look for a solution of the stochasticSchr¨odinger equation (5) of the form ψ t = ϕ t (cid:107) ϕ t (cid:107) (13)As in the completely positive divisible case, the change ofvariables (13) maps (5) into the linear problem d ϕ t = − ı H ϕ t − (cid:76) (cid:88) (cid:96) =1 (cid:18) Γ (cid:96),t † (cid:96) L (cid:96) d t − d ν (cid:96),t (L (cid:96) − (cid:19) ϕ t Once we know ϕ t , we can use (13) to determine the state vec-tor and the influence martingale. In particular, the influencemartingale is always factorizable into the form µ t = exp (cid:32) − (cid:90) t d s (cid:76) (cid:88) (cid:96) =1 (cid:103) (cid:96),t (cid:16) (cid:114) (cid:96),t (cid:107) L (cid:96) ψ t (cid:107) d t (cid:17)(cid:33) ˜ µ t (14)where ˜ µ t is a pure jump process. The factorization (14) maybe of use in numerical implementations. The description ofthe influence martingale method is thus complete. Interpretation.
The influence martingale process (10) ispositive definite whenever (cid:103) (cid:96),t ≥ − ∀ t and (cid:96) = 1 , . . . (cid:76) (15)These inequalities are well known in the theory of Poissonprocesses [43]. When they hold true, (9) simply correspondsto a redefinition of the rates of the Poisson processes. Inother words, at any time t the process µ t specifies the Radon–Nikodym derivative between two absolutely continuous Pois-son measures [25]. By (12), the violation of the inequalities(15) only occurs when (3) does not hold. In that case negativevalues of the influence process bring about a relative phase ( or π ) between instantaneous values of distinct quantum trajec-tories. We emphasize that the value of the phase is determinedonly by the individual trajectory and by the deterministic func-tions Γ (cid:96),t ’s encoding the state of the environment. No memoryeffects need to be taken into account. This is intuitively pleas-ant because in any finite dimensional Hilbert space (2) is justa matrix valued time non-autonomous linear ordinary differ-ential equation. Finally, it seems to us that the path integralformulation of quantum mechanics [44] offers a suggestiveaccount for the need to introduce the influence martingale.The influence martingale is a ”shadow on the cave wall” ofthe interference between paths in the Feynman path integraldescribing to the exact system environment unitary dynamics. ab FIG. 1:
Evolution of (cid:107) ϕ t (cid:107) (crosses, blue online) µ t (cid:107) ϕ t (cid:107) (circles, red online) and µ t (diamonds, green online). The initialdata ϕ = e where e is the H = σ + σ − . In the insets, Lamb shift(a) and Lindblad weight function (b) from [29] with β = − δ = 1 .The shaded regions in (b) emphasize negative values of Γ t . Example: Photonic band gap
The master equation inDirac’s interaction picture of a two level atom in a photonicband gap [29] is ˙ ρ t = S t ı [ σ + σ − , ρ t ] + Γ t ([ σ − ρ t , σ + ] + [ σ − , ρ t σ + ]) (16)where σ ± = ( σ ± ı σ ) / and { σ i } i =1 are Pauli matrices.The time dependent functions S t and Γ t are respectively theLamb shift and the Lindblad weight factor. Negative valuesof Γ t also imply a violation of the Kossakowski conditions([45] see also e.g, [17, 19]) a weaker form of positive divisi-bility than completely positive divisibility [46]. This fact ren-ders the unravelling of (16) in quantum trajectories particu-larly probing. In order to explore a genuine strong system-environment coupling, we proceed as in [32] and determine S t and Γ t by inserting equation (2.21) of [29] in the exact for-mulae (10.22), (10.23) of [2]. We also translate the origin oftime to t ≈ . so that S t and Γ t vanish at time origin. As [32]we perform the numerics with detuning δ = − β . In Fig. 1 weshow a typical realization of the squared norm of the process { ϕ t } t ≥ solution of the associated linear problem and howit is affected by the influence martingale. By (14) and in theabsence of jumps, the influence martingale amplifies exponen-tially the contributions of time intervals when Γ t is negative.When Γ t is positive µ t is constant in between jumps. A fur-ther, purely numerical challenge, presented by the unravellingof (16) stems from the fact that the ground state g of σ + σ − is preserved by the dynamics, whereas the Lindblad operatorcoincides with the annihilation operator. The consequence isa certain stiffness of the stochastic dynamics which is particu-larly evident for an initial datum ψ = e , e being the excitedstate. Fig. 2 reports the result of our numerical integrationsfor distinct values of the initial data. In all cases, influencemartingale weighed Monte Carlo averages over O (10 ) tra-jectories already well reproduce the non-monotonic evolutionof population probabilities within fluctuations induced errors. abcd FIG. 2:
Monte Carlo averages versus master equation predictions.The scatter plots are e † ρ t e for ψ = e (curve a , circles), and ψ = ( e + g ) / √ (curve c, diamonds) Curve b and d are ( e + g ) † ρ t ( e + g ) / for ψ = ( e + g ) / √ and ψ = ( e − g ) / √ , respectively. The continuous lines are thepredictions from the master equation. The shaded area of the size oftwice the square root of the variance of ensemble averages. Theensembles consists of trajectories for curves b,c, d. In the caseof the stiffer case of curve a, we verified that increasing theensemble up to × slowly damps the size of the fluctuations.The starting time mesh for the adaptive code is d t = 0 . . The occurrence ”quantum revivals” can be also quantitativelysubstantiated observing that the degree of memory effects in-troduced in [47] is simply related to Γ t for this model [17, 47].We performed the numerical integrations using the Tsitouras / Runge-Kutta method automatically switching for stiff-ness detection to to a -th order A-stable Rosenbrock. TheJulia code is offered ready for use in the ” DifferentialEqua-tion.jl ” open source suite [48]. Upon using [48], the numer-ical computations of Fig. 2 take a few minutes on a MBA11early 2014, definitely not a performance machine. Furthereasily conceivable numerical refinements (dedicated adaptiveschemes, larger than O (10 ) size of the ensemble, use of therepresentation (14) for increased efficiency) are beyond thepresent scopes. Example: A Λ -Qutrit is a three state { v i } i =1 system withenergies ω > ω > ω . In Dirac’s interaction picture thestate operator of Λ -qutrit satisfies (2) ( see e.g. [33]) with H = S ,t + S ,t v v † (17a) L (cid:96) = v (cid:96) v † , (cid:96) = 1 , (17b)We evaluate the Lamb shift S t and the Lindblad weights Γ (cid:96),t , (cid:96) = 1 , using the leading order time convolutionless pertur-bation theory expressions [2]: S (cid:96),t = γ λ ∆ (cid:96) (1 − e − λ t cos(∆ (cid:96) t )) − λ e − λ t sin(∆ (cid:96) t )∆ (cid:96) + λ Γ (cid:96),t = γ λ λ (1 − e − λ t cos(∆ (cid:96) t )) + ∆ (cid:96) e − λ t sin(∆ (cid:96) t )∆ (cid:96) + λ The numerical results of Fig. 3 exhibit a very good agreementbetween master equation and ensemble average predictions. . t v † ρ t v v † ρ t v v † ρ t v . ,t Γ ,t FIG. 3:
Population of the qutrit energy level associated to the initialcondition ψ = ( v + v ) / √ , with v i the eigenstate of the i -thenergy level. The black lines are the predictions of the masterequation specified by (17). The shaded area is determined by squareroot of the variance around the mean value of the influencemartingale weighed Monte Carlo average of N = 10 trajectoriesgenerated by (13). The left inset shows a zoomed in picture of theinitial population evolution. The right inset shows the rates withparameters γ = 1 . , λ = 0 . , ∆ = 8 λ and ∆ = 10 λ . The timemesh is d t = 0 . . Conclusions.
We have shown that master equations asso-ciated to divisible dynamical maps admit an unravelling intointerfering quantum trajectories. The influence martingale(10) amplifies or occasions cancellations between trajectories.As pointed out in [49], quantum trajectories can be inter-preted in three possible ways. The first is as mathematical tool to compute the solution of Lindblad–Gorini–Kossakowski–Sudarshan equation for high dimensional systems (see e.g.[10]). Second quantum trajectories are subjectively real : theirexistence and features are determined only contextually to agiven physical setup (see e.g. [23, 50]). And finally, theymight be an element of a still missing theory of quantum statereduction [12–15]. We believe that the interference propertieswe identify yield a significant contribution to, at least, the firsttwo of the interpretations mentioned above.
Acknowledgements
B.D. acknowledges support from theAtMath collaboration at the University of Helsinki. ∗ brecht.donvil@helsinki.fi[1] S. Haroche and J.-M. Raimond, Exploring the Quantum:Atoms, Cavities, and Photons (Oxford Graduate Texts, 2006)pp. X,616.[2] H.-P. Breuer and F. Petruccione,
The Theory of Open QuantumSystems , reprint ed. (Oxford University Press, 2002) pp. XXII,636. [3] H. M. Wiseman and G. J. Milburn,
Quantum Measurement andControl , 1st ed. (Cambridge University Press, 2009) p. 478.[4] T. Sauter, W. Neuhauser, R. Blatt, and P. E. Toschek, PhysicalReview Letters , 1696 (1986).[5] J. C. Bergquist, R. G. Hulet, W. M. Itano, and D. J. Wineland,Physical Review Letters , 1699 (1986).[6] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Rein-hold, R. Guti´errez-J´auregui, R. J. Schoelkopf, M. Mirrahimi,H. J. Carmichael, and M. H. Devoret, Nature , 1 (2019),arXiv:1803.00545 [quant-ph].[7] R. L. Hudson and K. R. Parthasarathy, Communications inMathematical Physics , 301 (1984).[8] A. Barchielli and V. P. Belavkin, Journal of Physics A: Math-ematical and General , 10.1088/0305-4470/24/7/022 (1991),quant-ph/0512189v1.[9] C. W. Gardiner, A. S. Parkins, and P. Zoller, Physical Review A , 4363 (1992).[10] J. Dalibard, Y. Castin, and K. M¨olmer, Physical Review Letters , 580 (1992).[11] H. Carmichael, An open systems approach to quantum optics:lectures presented at the Universit´e libre de Bruxelles, October28 to November 4, 1991 , Lecture Notes in Physics (Springer,1993).[12] N. Gisin, Physical Review Letters , 1657 (1984).[13] G.-C. Ghirardi, A. Rimini, and T. Weber, Physical Review D , 470 (1986).[14] A. Bassi and G. Ghirardi, Physics Reports , 257–426(2003).[15] I. C. Percival, Quantum State Diffusion (Cambridge UniversityPress, 2003).[16] K. Kraus,
States, Effects, and Operations Fundamental Notionsof Quantum Theory , edited by K. Kraus, A. B¨ohm, J. D. Dol-lard, and W. H. Wootters (Springer Berlin Heidelberg, 1983).[17] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Reviews ofModern Physics , 021002 (2016), arXiv:1505.01385 [quant-ph].[18] E. B. Davies, Quantum Theory of Open Systems (AcademicPress, 1976) p. 171.[19] A. Rivas and S. F. Huelga,
Open Quantum Systems ,Springer Briefs in Physics (Springer Berlin Heidelberg, 2012)arXiv:1104.5242 [quant-ph].[20] G. Lindblad, Communications in Mathematical Physics , 119(1976).[21] V. Gorini, A. Kossakowski, , and E. C. G. Sudarshan, Journalof Mathematical Physics , 821 (1976).[22] D. Chru´sci´nski and S. Pascazio, Open Systems & InformationDynamics , 1740001 (2017), arXiv:1710.05993 [quant-ph].[23] H.-P. Breuer and F. Petruccione, Physical Review E ,428–441 (1995).[24] A. Barchielli and C. Pellegrini, Journal of MathematicalPhysics 10.1063/1.3514539 (2010), 1006.4527v1.[25] F. C. Klebaner, Introduction to stochastic calculus with appli-cations , 2nd ed. (Imperial College Press, 2005) p. 416. [26] R. P. Feynman and F. L. J. Vernon, Annals of Physics , 118(1963).[27] M. W. Y. Tu and W.-M. Zhang, Physical Review B , 235311(2008), arXiv:0809.3490 [cond-mat.mes-hall].[28] B. Donvil, D. Golubev, and P. Muratore-Ginanneschi, PhysicalReview B , 245401 (2020), 2007.15923.[29] S. John and T. Quang, Physical Review A , 1764 (1994).[30] H.-P. Breuer, B. Kappler, and F. Petruccione, Physical ReviewA , 1633 (1999), arXiv:quant-ph/9906024 [quant-ph].[31] H.-P. Breuer, Physical Review A , 012106 (2004), quant-ph/0403117.[32] J. Piilo, S. Maniscalco, K. H¨ark¨onen, and K.-A. Suominen,Physical Review Letters , 180402 (2008), arXiv:0706.4438[quant-ph].[33] J. Piilo, S. Maniscalco, K. H¨ark¨onen, and K.-A. Suominen,Physical Review A , 10.1103/physreva.79.062112 (2009),0902.3609.[34] A. Smirne, M. Caiaffa, and J. Piilo, Physical Review Letters10.1103/physrevlett.124.190402 (2020), 2004.09537.[35] L. Di´osi and W. T. Strunz, Physics Letters A , 569 (1997).[36] L. Di´osi, N. Gisin, and W. T. Strunz, Physical Review A ,1699 (1998), arXiv:quant-ph/9803062 [quant-ph].[37] J. Gambetta and H. M. Wiseman, Phys. Rev. A ,10.1103/physreva.66.012108 (2002).[38] A. Tilloy, Quantum , 29 (2017), 1701.01948v3.[39] M. Chen and J. Q. You, Physical Review A , 052108 (2013),arXiv:1203.2217 [quant-ph].[40] J. Gambetta and H. M. Wiseman, Physical Review A ,062104 (2003), quant-ph/0307078.[41] J. Gambetta and H. M. Wiseman, Foundations of Physics ,419 (2004).[42] M. W. Jack and M. J. Collett, Physical Review A ,10.1103/physreva.61.062106 (2000).[43] R. K. Boel, P. Varaiya, and E. Wong, SIAM Journal on Control , 1022 (1975).[44] R. P. Feynman and A. R. Hibbs, Quantum mechanics and pathintegrals: Emended Edition , edited by D. F. Styer (Dover Pub-lications, 2010) p. 384.[45] A. Kossakowski, Reports on Mathematical Physics , 247(1972).[46] D. Chru´sci´nski and S. Maniscalco, Physical Review Letters , 10.1103/physrevlett.112.120404 (2014).[47] H.-P. Breuer, E.-M. Laine, and J. Piilo, Physical Review Letters , 210401 (2009), arXiv:0908.0238 [quant-ph].[48] C. Rackauckas and Q. Nie, Journal of Open Research Software (2017).[49] H. M. Wiseman, Quantum and Semiclassical Optics: Journal ofthe European Optical Society Part B , 205 (1996).[50] A. H. Wilson and G. J. Milburn, Physical Review A47