Stochastic unravellings of non-Markovian completely positive and trace preserving maps
aa r X i v : . [ qu a n t - ph ] O c t Stochastic unravellings of non-Markovian completely positiveand trace preserving maps
G. Gasbarri ∗ and L. Ferialdi
2, 3, † Department of Physics, University of Trieste, Strada costiera 11, 34151 Trieste, Italy Department of Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia Department of Physics and Astronomy, University of Southampton, SO17 1BJ, United Kingdom (Dated: October 11, 2018)We consider open quantum systems with factorized initial states, providing the structure of thereduced system dynamics, in terms of environment cumulants. We show that such completelypositive (CP) and trace preserving (TP) maps can be unraveled by linear stochastic Schrödingerequations (SSEs) characterized by sets of colored stochastic processes (with n -th order cumulants).We obtain both the conditions such that the SSEs provide CPTP dynamics, and those for unravelingan open system dynamics. We then focus on Gaussian non-Markovian unravellings, whose knownstructure displays a functional derivative. We provide a novel description that replaces the functionalderivative with a recursive operatorial structure. Moreover, for the family of quadratic bosonicHamiltonians, we are able to provide an explicit operatorial dependence for the unravelling. I. INTRODUCTION
Stochastic Schrödinger equations (SSEs) were intro-duced in the framework of measurement theory, wherethe effect of repeated measurements is described by theintroduction of a Wiener process in the Schrödinger equa-tion [1]. In the following years, different SSEs were pro-posed in the context of collapse models where one mod-ifies the Schrödinger dynamics in such a way to accountfor the wave packet reduction within a unique dynam-ical principle [2]. Markovian SSEs, i.e. those display-ing white noises, have been deeply investigated thanks tothe possibility of exploiting Ito stochastic calculus [3–5].The solutions of these SSEs “unravel” (i.e. reproduce inaverage) Markovian open systems dynamics. However,many physical systems cannot be described by Marko-vian dynamics, and they need non-Markovian ones [6, 7].This fact pushed the development and investigation ofnon-Markovian SSEs. The extension to non-Markovianquantum state diffusion was first tackled by Diosi andStrunz in [8], where Gaussian, complex, colored stochas-tic processes were considered. These seminal works werefollowed by others both in the realm of collapse mod-els [9–12], and in that of continuous measurement [13](for a more detailed historical review see [7]). We shouldmention that in this context one considers stochasticequations with discontinuous or jump processes, both forMarkovian and non-Markovian dynamics [2, 14]. We areinstead interested in stochastic processes that are con-tinuous and continuously differentiable, and we will fo-cus on these in the remainder of the paper. SSEs arealso exploited, as method alternative to master equa-tions, to numerically investigate open quantum systems.The advantage is that SSEs typically require less numer- ∗ Electronic address: [email protected] † Electronic address: [email protected] ical resources than master equations (because the wavefunctions Hilbert spaces have smaller dimensionality thanthose of density matrices). Different expansions methodshave been proposed to investigate non-Markovian SSEs,but only in few cases these equations have been analyti-cally solved [10, 11, 15, 16]. We should also mention theoriginal method proposed by Tilloy [17], who introducesnew stochastic terms in order to rewrite non-MarkovianSSEs as averages over explicit time-local (bi-)stochasticdifferential equations. Recently, the first full character-ization of Gaussian linear stochastic unravellings havebeen provided in [18]. Such a characterization is im-plicit, as it depends on a functional derivative (furtherexplained in Sec.IV). The symmetries of these unravel-lings were analyzed in [19], while a measurement inter-pretation in terms of Bargmann states has been providedin [20].Most of the results on non-Markovian SSEs in the lit-erature concern Gaussian stochastic processes, and onlyfew papers attempt to use continuous, non-Gaussianstochastic processes in very specific cases [21]. The aimof this paper is to take a step forward in the theoreticaldescription of non-Markovian unravelings. On the oneside, we show that it is possible to unravel non-Gaussian,completely positive (CP) and trace preserving (TP) non-Markovian open system dynamics, and we provide thecondition for the existence of such an unraveling (Sec.III).In general, if the open system dynamics is non-Gaussian,the SSE will display continuous colored stochastic pro-cesses characterized by their n-th order cumulants. Onthe other hand, we provide a novel description of Gaus-sian non-Markovian unravelings that does not make useof the functional derivative and explicitly depends on thesystem operators (Sec.IV). In the next section we start bywe investigating the structure of (non-Gaussian) CPTPopen system dynamics, that is the reference dynamicsthat we want to unravel.
II. OPEN QUANTUM SYSTEMS MAP
We consider a system ( S ) interacting with a genericenvironment ( E ). The evolution of the open system den-sity matrix ˆ ρ SE in the interaction picture is described bythe Von-Neumann equation ( ~ = 1 ) i∂ t ˆ ρ SE ( t ) = h ˆ V t , ˆ ρ SE ( t ) i , (1)where ˆ V t is a generic interaction Hamiltonian betweenthe system and the environment, that can be rewrittenin total generality as: ˆ V t = X α ˆ f α ( t ) ˆ φ α ( t ) , (2)where ˆ f α ( t ) and ˆ φ α ( t ) respectively are Hermitian systemand environment operators. The solution of Eq. (1) canbe formally written as: ρ SE ( t ) = ˆΦ t ρ SE ˆΦ † t (3)with ˆΦ t = T exp " − i Z t dτ X α ˆ f α ( τ ) ˆ φ α ( τ ) . (4)The time ordering T acts on a generic function of thetype f ( R t dτ ˆ V τ ) , by ordering the products of operators ˆ V τ of its Taylor expansion: T f (cid:18)Z t dτ ˆ V τ (cid:19) = ∞ X n =0 n ! T "(cid:18)Z t dτ ˆ V τ (cid:19) n ∂ n f ( x ) ∂x n (cid:12)(cid:12)(cid:12)(cid:12) x =0 , (5)provided the series exists. Since we are interested in theeffective evolution of the system S we trace over the en-vironment degrees of freedom to obtain: ˆ ρ S ( t ) = Tr E [ ˆΦ t ˆ ρ SE ˆΦ † t ] . (6)Under the assumption of factorized initial state ˆ ρ SE =ˆ ρ S ⊗ ˆ ρ E , the effective evolution can be described by theaction of a dynamical map f M t on the system initial state ρ s . By replacing Eq. (4) in Eq. (6), one can write thedynamical map as: f M t [ˆ ρ S ] = D T e − i R t dτ P α [ f + α ( τ ) φ − α ( τ )+ f − α ( τ ) φ + α ( τ )] ˆ ρ S E (7)where h . . . i = Tr E [ . . . ˆ ρ E ] , and the superscripts + / − de-note the super-operators f ± ˆ ρ = ˆ f ˆ ρ ± ˆ ρ ˆ f (similar defini-tions hold for the super-operators φ ± ). We rewrite themap f M t as a time ordered exponential with a time de-pendent generator, by means of a cumulant expansion(see Appendix A for explicit calculation): f M t [ˆ ρ S ] = T exp " ∞ X n =1 ( − i ) n Z t d τ n ˜ k n ( τ n ) ˆ ρ S , (8) where τ n = ( τ , · · · τ n ) ( n denoting the length of thearray). The terms of the expansion are defined as ˜ k n ( τ n ) = n X q =0 X α i , P q f − α ( τ ) f l α ( τ ) ··· f l n α n ( τ n ) e C +¯ l ··· ¯ l n α ··· α n ( τ n ) , (9)where e C are the ordered bath cumulants defined inEq. (A9). We used the following notation: l i ∈ { + , −} , ¯ l i = − l i , P q is the permutation of the indexes l i such that q is the number of plus signs in the array ( l , · · · , l n ) , and α i denotes the sum over all { α , · · · , α n } . We stress thatthe first super-operator on the left in ˜ k n ( τ n ) is always f − .This is an important feature because it guarantees thatthe map is trace preserving (indeed Tr S [ f − ˆ O ˆ ρ s ] = 0 ).Note moreover that the time ordering acts on the ba-sic elements f , hence one has first to replace Eq. (9) inEq. (8) , and then apply T .We eventually notice that Eq. (8) recovers the one de-rived by Diosi and Ferialdi [18] when the (Gaussian) en-vironment is completely characterized by its second cu-mulant. Indeed, in such a case Eq. (8) reduces to f M t = T exp " − Z t d τ X α ,α f − α ( τ ) × (10) × (cid:16) f − α ( τ ) e C ++ α α ( τ ) + f + e C + − α α ( τ ) (cid:17) , that coincides with Eq.(17) in [18]. III. GENERAL UNRAVELLING
We now consider a stochastic Schrödinger equation ininteraction picture i∂ t | ψ t i = ˆ V t | ψ t i . (11)where ˆ V t is an operator valued stochastic field in theHilbert space H , that can be decomposed as a sum ofHermitian system operators ˆ f α like in Eq. (2), where theoperators ˆ φ α are replaced by a set of complex stochas-tic processes { φ α } . A decomposition with non-Hermitianoperators can be obtained performing a unitary transfor-mation as shown in [19]. Unlike previous works on non-Markovian stochastic unravellings, we consider generalstochastic processes, that are completely characterizedby their n -th order correlation functions, with n arbi-trary. Equation (11) is solved by Eq. (4) with ˆ φ → φ ,and the average dynamics is described by the map M t ,defined as follows: M t [ˆ ρ S ] = D ˆΦ t ˆ ρ S ˆΦ † t E (12)where ˆ ρ S = | ψ i h ψ | is a pure initial state, and now h· · ·i denotes the stochastic average. In order for M t tobe a physical map, it must be CPTP. The definition (12)itself can be understood as the Kraus-Stinespring decom-position of M t , guaranteeing its CP [22]. By replacingEq. (4) in Eq. (12), one can conveniently rewrite M t likein Eq. (7), where now φ ± = φ ± φ ∗ are twice the realand imaginary parts of φ . In order to find the condi-tions under which the averaged map is TP, we perform–analogously to the quantum case– a cumulant expansion(see Appendix A), obtaining Eq. (8) where ˜ k is replacedby k n ( τ n ) = n X q =0 X α i , P q f l α ( τ ) ··· f l n α n ( τ n ) C ¯ l ··· ¯ l n α ··· α n ( τ n ) . (13) C ¯ l ··· ¯ l n α ··· α n ( τ n ) are the n -th order cumulants of the set ofstochastic processes φ α (see Appendix A for the explicitexpression). It is important to stress that the stochasticcumulants C in principle differ from the quantum cumu-lants e C . This is due to the fact that the trace over quan-tum degrees of freedom in general cannot be describedas a stochastic average.We perform the trace of M t , and we observe that thecontributions of the type f − α ( τ ) · · · f l n α n ( τ n ) C + ··· ¯ l n α ··· α n ( τ n ) (i.e. those with l = − ) in each k n ( τ n ) are killed by thecyclicity of the trace (remember that f − is a commutator,and accordingly Tr { f − ˆ O } = 0 ): Tr {M t [ˆ ρ S ] } =Tr T exp Z t d τ n ∞ X n =1 n X q =1 X α i , P q ( − i ) n f + α ( τ ) f l α ( τ ) ··· f l n α n ( τ n ) C − ¯ l ··· ¯ l n α ··· α n ( τ n ) ˆ ρ S ) . (14)Recalling that M t is TP if Tr {M t [ˆ ρ S ] } = Tr { ˆ ρ S } , onefinds that this condition is satisfied only if the exponentin the previous equation is zero, i.e. if either the series oreach of its terms are zero. In the first case one needs tofine-tune the stochastic cumulants C (by manipulatingthe processes φ α ) in such a way that the series in theexponent sum to zero. This is in general a very difficulttask and we are not aware of any technique that allowsto provide an explicit constraint. The second option,though being a more stringent requirement, allows to findand explicit constraint: C − ¯ l ··· ¯ l n α ··· α n ( τ n ) = 0 , ∀ n , (15)i.e. all the cumulants obtained tracing at least one purelyimaginary noise φ − are zero. Accordingly, the averagemap M t is TP if it is generated by an interaction po-tential ˆ V ( t ) that displays only real stochastic processes φ α . In this case ˆ V ( t ) is purely Hermitian, and Eq. (13)is replaced by ( φ α = φ + α ) k n ( τ n ) = X α i f − α ( τ ) ··· f − α n ( τ n ) C + ··· + α ··· α n ( τ n ) . (16) However, it is well known that a SSE of this kind gen-erates a dynamics that cannot describe dissipative phe-nomena [10]. Dissipative dynamics are indeed unraveledby SSEs displaying an interaction potential (2) with com-plex stochastic processes [11, 23] (provided that the sys-tem operators are Hermitian). Since the previous cal-culations explicitly show that the unravelling (11) withcomplex noises leads to a dynamics that is not TP, wenow investigate whether it is possible to obtain such a TPmap starting from a SSE different from (11). We do so byadding new terms to the stochastic map ˆΦ t , and we seekthe conditions on these new terms such that the aver-age map is TP. We stress that one might make M t tracepreserving by modifying directly the k n ( τ n ) of Eq. (13).This operation, though leading to a TP map, does notguarantee that M t preserves the factorized structure ofEq. (12), which is an essential requirement to obtain aSSE unravelling it.The most general modification of ˆΦ t is obtainedby adding to the exponent an operator functional g ( R t dτ ˆ f τ ) , that for later convenience we Taylor expandaround R t dτ ˆ f τ = 0 (see e.g. Appendix A of [24]). Thisleads to a new map ˆΞ t : ˆΞ t = T exp " − i Z t dτ X α ˆ f α ( τ ) φ α ( τ ) − ∞ X n =1 Z t d τ n ˆ O n ( τ n ) , (17)where the operators ˆ O n ( τ n ) are defined by Z t d τ n ˆ O n ( τ n ) = 1 n ! (cid:18)Z t dτ ˆ f τ (cid:19) n ∂ n g ( x ) ∂x n (cid:12)(cid:12)(cid:12)(cid:12) x =0 . (18)We stress that ˆ O n ( τ n ) are deterministic: if they werestochastic they would lead to a new map that could berewritten in the form (4) (with new stochastic processes),therefore not solving the trace issue. The stochas-tic map (17) generates the average dynamics N t [ˆ ρ S ] = h ˆΞ t ˆ ρ S ˆΞ † t i . We remark that this dynamics is CP evenif ˆΞ t is truncated at the n -th order, because the Kraus-Stinespring structure of Eq. (12) is preserved. Exploitingthe identity T e i ( a + b ) = T e i T log ( e ia )+ b and the calcula-tions in Appendix A, one finds that N t = T exp " ∞ X n =1 Z t d τ n (cid:16) ( − i ) n k n ( τ n ) − O n ( τ n ) (cid:17) (19)with O n ( τ n ) = H + n ( τ n ) + A − n ( τ n ) . (20)The super-operators H + n , A − n are built –with the knownrules– respectively from the Hemitian and anti-Hermitianparts of ˆ O n : ˆ H n = ( ˆ O n + ˆ O † n ) / and ˆ A n = ( ˆ O n − ˆ O † n ) / .By performing the trace of Eq. (19), one finds that theaveraged map N t is trace preserving if the following con-dition is satisfied (see Appendix B): ˆ A n ( τ n ) =( − i ) n n X q =1 X α i , P q ˆ f α ( τ ) f ↼ l α ( τ ) · · · f ↼ l n α n ( τ n ) C − ¯ l ··· ¯ l n α ··· α n ( τ n ) , (21)where the super-operators f ↼ ± are defined as ˆ O f ↼ ± =ˆ O ˆ f ± ˆ f ˆ O , the arrow denoting the fact that these super-operators are acting on the left. We recall that ¯ l = − because of the trace operation (see Eq. (14)). We thenfind that the request of N t being TP map can be satis-fied –for a generic noise– if a multi-time anti-Hermitianoperator of the form in Eq. (21) is added to the map ˆΞ t .Accordingly, by replacing Eq. (21) in Eq. (17) we obtainthe structure of a general stochastic map that generatesan average CPTP dynamics. Interestingly, the fact thatonly the anti-Hermitian ˆ A n are involved in the TP condi-tion implies that the Hermitian contributions ˆ H n can bechosen freely. This degree of freedom can be exploited totune the average dynamics N t obtained from a stochasticevolution, in such a way that it recovers an open systemdynamics f M t . By matching the exponents of Eq. (19)and Eq. (8), one finds that ˆΞ t unravels f M t if H + n ( τ n ) = ( − i ) n (cid:16) k n ( τ n ) − ˜ k n ( τ n ) (cid:17) − A − n ( τ n ) . (22)One might question whether it is possible to find a moreexplicit expression for this formal equation, e.g. by de-composing the operator ˆ O n over the set { ˆ f α } in Eq. (17).Unfortunately, such a decomposition is helpful only in theGaussian case. The calculations showing this are ratherinstructive and we reported them in Appendix B. As amatter of facts, there are two cases for which one can ob-tain a simpler expression for ˆ H n : when all operators ˆ f α ( t ) commute at any time ( [ ˆ f α ( t ) , ˆ f α ( τ )] = 0 ), and the Gaus-sian case (stochastic processes characterized only by theirfirst two cumulants). The first case is trivially solved andwe will not report it here, while the Gaussian case is moreinteresting and it will be considered it in the next section. IV. GAUSSIAN UNRAVELLING
We consider complex Gaussian stochastic processes,i.e. those completely characterized by their first two cu-mulants: h φ ± α ( t ) i = 2 C ± α ( t ) h φ ± α ( t ) ∗ φ ± β ( τ ) i − (cid:10) φ ± α ( t ) ∗ (cid:11) h φ ± β ( τ ) i = 4 C ±± αβ ( t, τ ) . (23) The CPTP stochastic map (17) associated to such pro-cesses simplifies to ˆΞ t = T exp (cid:26) − i Z t dτ h ˆ f α ( τ ) (cid:2) φ α ( τ ) − C − α ( τ ) (cid:3) + Z t d τ ˆ f α ( τ ) h f ↼ β + ( τ ) C −− αβ ( τ )+ f ↼ β − ( τ ) C − + αβ ( τ ) i − Z t dτ ˆ H ( τ ) − Z t d τ ˆ H ( τ ) (cid:27) , (24)where we are now making use of the Einstein conven-tion of summing over repeated indexes. In this specificcase we can set ˆ H ( τ ) = 0 , and determine the operator ˆ H ( τ ) in such a way that the average map reproducesthe effective Gaussian map generated by the trace overan environment: ˆ H ( τ ) = ˆ f α ( τ ) h f ↼ − β ( τ ) C −− αβ ( τ ) + f ↼ + β ( τ ) C − + αβ ( τ ) i . (25)Replacing this equation in Eq. (26) we obtain the follow-ing expression for the stochastic map ˆΞ t = T exp (cid:26) − i Z t dτ ˆ f α ( τ ) (cid:0) φ α ( τ ) − C − α ( τ ) (cid:1) (26) + Z t d τ ˆ f α ( τ ) ˆ f β ( τ ) (cid:16) C −− αβ ( τ ) + C − + αβ ( τ ) (cid:17)(cid:27) . Since C − α ( τ ) in (26) simply gives a shift of the mean valueof φ α , we absorb it as follows φ α ( τ ) − C − α ( τ ) → φ α ( τ ) ,in order to simplify the notation. This leads to ˆΞ t = T n e − i R t dτ ˆ f α ( τ ) [ φ α ( τ ) − i R t dτ ˆ f β ( τ ) K αβ ( τ,τ ) ] o , (27)with the following conditions on the new φ : h φ α ( t ) i = h φ α ( t ) ∗ i K αβ ( t, τ ) = 12 (cid:0) h φ α ( t ) ∗ φ β ( τ ) i − h φ α ( t ) φ β ( τ ) i (cid:1) θ tτ . (28)Equation (27) is the most general Gaussian, non-Markovian stochastic map that unravels a CPTP dissipa-tive dynamics, and coincides with the one first obtainedin [18] for h φ α i = 0 . In order to find the SSE associatedto this map, we need to differentiate it with respect to t ,obtaining ∂ t ˆΞ t = − i ˆ f α ( t ) φ α ( t ) ˆΞ t (29) − T (cid:20) ˆ f α ( t ) Z t dτ K αβ ( t, τ ) ˆ f β ( τ ) × e − i R t dτ ˆ f α ( τ ) [ φ α ( τ ) − i R τ dτ ˆ f β ( τ ) K αβ ( τ,τ ) ] i . The main issue with this equation is that the operator ˆ f β ( τ ) in the second line is entangled with the exponentialthrough the time ordering. The approach most widelyused to deal with this issue is to exploit the Furutsu-Novikov theorem [25], that allows to rewrite Eq. (29) asfollows: ∂ t ˆΞ t = − i (cid:20) ˆ f α ( t ) φ α ( t ) − i ˆ f α ( t ) Z t dτ K αβ ( t, τ ) δδφ β ( τ ) (cid:21) ˆΞ t , (30)where δ/δφ β ( τ ) denotes a functional derivative. How-ever, this is just a formal and elegant rewriting ofEq. (29), that it is very difficult to exploit for practi-cal purposes [10, 11, 15, 26]. We propose an alternativeperturbative scheme, that has the merit of allowing for arecursive definition of the expansion terms for the timelocal generator. We rewrite Eq. (29) as follows ∂ t ˆΞ t = − i h ˆ f α ( t ) φ α ( t ) + ˆΛ t i ˆΞ t (31)where ˆΛ t is a time-local representation of the non-Markovian contribution to the dynamics. When inverted,this equation allows to write a formal expression for ˆΛ t : ˆΛ t = h i ∂ t ˆΞ t − ˆ f α ( t ) φ α ( t ) ˆΞ t i ˆΞ − t . (32)We adopt the same strategy as in [27]: we Taylor expandthe stochastic map (26) to obtain ˆΞ t = ∞ X n =0 n X k =0 n ! k ! T (cid:26) (cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n − k × (cid:18) − Z t d τ ˆ f α ( τ ) ˆ f β ( τ ) K αβ ( τ ) (cid:19) k (cid:27) , (33)then after some manipulation and explicitly solving thetime ordering (see Appendix D), we write the map ˆΞ t asfollows: ˆΞ t = ∞ X n =0 ( − i ) n ˆ ξ n (34)with ˆ ξ n = Z t d τ n X α i ˆ f α ( τ ) ··· ˆ f α n ( τ n ) × ⌊ n/ ⌋ X k =1 X P ( a ) K α α ( τ a , τ a ) ··· K α k − α k ( τ a k − , τ a k ) × φ α k +1 ( τ a k +1 ) ··· φ α n ( τ a n ) θ τ a k +1 ··· τ an , (35)where ⌊ n/ ⌋ denotes the floor function of n/ (i.e. thegreatest integer less than or equal to n/ ), and P ( a ) isthe permutation of the indexes a i ∈ { , . . . , n } . We alsoexpand the right hand side of Eq. (32) (see Appendix Dfor detailed calculation.): h i ∂ t ˆΞ t − ˆ f α ( t ) φ α ( t ) ˆΞ t i = ∞ X n =1 ( − i ) n ˆ d n +1 (36) where ˆ d n +1 = X α i ˆ f α ( t ) Z t d τ n ˆ f α ( τ ) ··· ˆ f α n ( τ n ) × ⌊ n/ ⌋ X k =1 X P ( a ) K αα ( t, τ a ) ··· K α k − α k − ( τ a k − , τ a k − ) | {z } k elements × φ α k ( τ a k ) ··· φ α n ( τ a n ) | {z } ( n − k ) elements θ τ a k ··· τ an . (37)Exploiting the result in Appendix B of [27] we can nowrewrite the inverse map as: ˆΞ − t = ∞ X n =0 ( − i ) n ˆ M n , (38)with M = 1 and ˆ M n = n X k =1 ˆ M k ˆ ξ n − k . (39)Recollecting these results, replacing them in Eq. (32),and following the strategy in Appendix C of [27], weeventually obtain the following perturbative series for thegenerator ˆΛ t : ˆΛ t = ∞ X n =1 ˆ L n (40)whose elements ˆ L n are defined by the following recursiveformula ˆ L n = ˆ d n − n X k =1 ˆ L k ˆ ξ n − k . (41)It is important to stress that the order n in the expan-sion terms ˆ ξ n , ˆ d n and ˆ L n denotes the power of operators ˆ f α displayed by them. This allows to understand theseries as an expansion in the coupling strength of the in-teraction, which provides a useful tool for the practicalanalysis of non-Markovian unravellings. Equation (30)instead, thought being much more elegant cannot be di-rectly used –unless the functional derivative is known–for explicit calculations. The explicit form of the ˆ L n isstill rather involved, making the evaluation of higher or-ders demanding. However, for a specific class of physicalsystems this problem can be solved and an explicit un-ravelling obtained. A. Bosonic quadratic Hamiltonian
We assume that the system of interest is bosonic and itsfree dynamics is described by a quadratic Hamiltonian.The advantage provided by this family of systems is thatthey allow to apply Wick’s theorem to disclose the timeordering of Eq. (29). Indeed, these systems display lin-ear Heisenberg equations of motions, and accordingly thecommutators [ ˆ f α ( τ ) , ˆ f β ( s )] are c-functions. An exampleof system falling in this category is a damped harmonicoscillator: in this case the operators ˆ f are simply positionand momentum of the oscillators, and contractions arecombination of sine and cosine functions [28]. The strat-egy we adopt is the following: we start from Eq. (33)and, instead of solving explicitly the time ordering, wekeep its formal expression. Differentiating Eq. (33) andexploiting the properties of the Cauchy product of twoseries and we obtain: i∂ t ˆΞ n = ˆ f j ( t ) φ j ( t ) ˆΞ t (42) − i ˆ f α ( t ) T "Z t dτ K αβ ( t, τ ) ˆ f β ( τ ) ∞ X n =0 ∞ X k =0 · · · where the dots denote the same argument of the series ofEq. (33). Note that this is just a convenient rephrasingof Eq. (29). Our aim is to achieve an equation of thetype (31), i.e. to express ∂ t ˆΞ t in terms of ˆΞ t . By applyingWick’s theorem we find that the time ordering in thesecond line of Eq. (42) can be rewritten as follows: T (Z t dτ K αβ ( t, τ ) ˆ f β ( τ ) ∞ X n =0 ∞ X k =0 ··· ) = Z t dτ K αβ ( t, τ ) (cid:18) ˆ f β ( τ ) − i Z t dτ C βα ( τ, τ ) φ α ( τ ) (cid:19) Ξ t − T (Z t dτ ˆ f β ( τ ) K (2) αβ ( t, τ ) ∞ X n =0 ∞ X k =0 · · · ) (43)where K (2) is a suitably defined kernel, and C βα ( τ , τ ) =[ ˆ f β ( τ ) , ˆ f α ( τ )] θ τ τ is a Wick contraction (see AppendixE). We have been able to decompose the initial time or-dering in two terms: one proportional to ˆΞ t , and one dis-playing a time ordering that has the same structure asthe initial one, but with a different kernel. We can thenapply Eq. (43) to this new time ordering, and iteratingthis procedure we obtain (see Appendix E) ∂ t ˆΞ t = (cid:20) − i ˆ f α ( t ) φ α ( t ) − i ˆ f α ( t ) Z t dτ K αβ ( t, τ ) (44) × (cid:18) − i ˆ f β ( τ ) + Z t dτ C βα ( τ, τ ) φ α ( τ ) (cid:19)(cid:21) ˆΞ t , where K = P ∞ n =1 K ( n ) , with K (1) = K , and K ( n ) definedrecursively through the following formula K ( n ) αβ ( t, τ ) = Z t d τ K ( n − αα ( t, τ ) C α α ( τ , τ ) K α β ( τ , τ ) , (45) with K α β ( τ , τ ) = K α β ( τ , τ ) + K βα ( τ, τ ) . The con-dition for convergence of the series is: lim n →∞ K ( n ) αβ ( t, τ ) → . (46)This equation provides the explicit non-Markovian un-ravelling for the family of bosonic quadratic systems, gen-eralizing the result obtained in [10] for a specific model.Eventually, the average dynamics associated to this equa-tion is described by the master equation obtained in [28]. V. CONCLUSIONS
We have first investigated open quantum systems withfactorized initial states, and we derived the structure ofthe reduced CPTP map, making explicit its dependenceon bath cumulants. We then considered the class of linearSSEs with continuous stochastic processes, exploring thepossibility of using them to construct CPTP dynamics.We found that when real stochastic process are consid-ered the structure of the unravelling is trivial. On theother hand, we showed that when SSEs with complexnoises are considered, the conditions such that these pro-vide average CPTP maps are rather involved. Further-more, we provided the condition that SSEs have to fulfilin order to unravel open quantum systems dynamics.We eventually focused on Gaussian unravellings, pro-viding a new perturbative expansion that relies on a re-cursive equation. This expansion, though being less ele-gant than the known one displaying a functional deriva-tive, allows to compute recursively the non-local term ofthe unravelling. Furthermore, for the family of quadraticbosonic systems, we have been able to provide an explicitoperatorial expression for the unravelling.
Acknowledgments
The authors are grateful to L. Diosi for spotting aflaw on a first derivation of the conditions for the non-Gaussian stochastic unraveling. The authors furtherthank A. Bassi for useful discussions. GG acknowl-edges funding from the European Union’s Horizon 2020research and innovation programme under grant agree-ment No 766900. LF acknowledges funding from theRoyal Society under the Newton International Fellow-ship No NF170345. During his stay at the University ofLjubljana, LF has been supported by the ERC AdvancedGrant 694544 - OMNES.
Appendix A
In this Appendix we provide explicit calculation leading to Eqs. (8)-(9). We start from Eq. (7): M t [ˆ ρ S ] = D T e − i R t dτ P α [ f + α ( τ ) φ − α ( τ )+ f − α ( τ ) φ + α ( τ )] ˆ ρ S E (A1)where f ± ( φ ± ) are super-operators acting respectively on the system (environment) Hilbert space (definition in themain text), and h· · ·i denotes the average over a certain measurable space (in this case the trace over the environmentdegrees of freedom). Exploiting Eq. (5) we can define the time ordered logarithm as follows: T log (cid:20) f (cid:18)Z t dτ ˆ V τ (cid:19)(cid:21) = ∞ X n =1 ( − ) n +1 n T " f (cid:18)Z t dτ ˆ V τ (cid:19) n (A2)provided that | f ( R t ˆ V τ ) | < . One can check that such time ordered logarithm satisfies the identity T log( T e R t dτ ˆ V τ ) = R t dτ ˆ V τ . Exploiting this in Eq. (A1) we find: M t = T exp (cid:26) T log D T e − i R t dτ P α [ f + α ( τ ) φ − α ( τ )+ f − α ( τ ) φ + α ( τ )] E(cid:27) . (A3)The exponent of this expression can be rewritten as the following limit lim ε ± α ( τ ) → T log D T e − i R t dτ [( f + α ( τ )+ iε + α ( τ )) φ − α ( τ )+( f − α ( τ )+ iε − α ( τ )) φ + α ( τ )] E , (A4)and exploiting functional calculus [30], one can conveniently rewrite the logarithm in Eq. (A3) as follows T log D T e − i R t dτ P α [ f + α ( τ ) φ − α ( τ )+ f − α ( τ ) φ + α ( τ )] E = lim ε ± α ( τ ) → T e − i R dτ P α [ f + α ( τ ) δ/δε − α ( τ )+ f − α ( τ ) δ/δε + α ( τ )] log D T e R t dτ P β [ ε − β ( τ ) φ − β ( τ )+ ε − β ( τ ) φ + β ( τ )] E (A5)where δ/δε ± α ( τ ) denotes a functional derivative. This expression is convenient because it allows to evaluate indepen-dently the system and the environment contributions to the system dynamics. One may notice indeed that, in ther.h.s. of Eq. (A5), the super-operators f ± acting on the system are all collected in the first exponential, while theremaining degrees of freedom are displayed only by the logarithm. Expanding the exponential function of the systemsuper-operators f ± α ( τ ) and resolving its time ordering one obtains T log D T e − i R t dτ P α [( f + α ( τ )+ iε − α ( τ )) φ − α ( τ )+( f − α ( τ )+ iε − α ( τ )) φ + α ( τ )] E = ∞ X i =0 ( − i ) n n ! Z t d τ n T (cid:26) X α (cid:20) f + α ( τ ) δδε − α ( τ ) + f − α ( τ ) δδε + α ( τ ) (cid:21) n (cid:27) log D T e R t dτ P β [ ε − β ( τ ) φ − β ( τ )+ ε − β ( τ ) φ + β ( τ )] E = ∞ X i =0 ( − i ) n Z t d τ n n X q =0 X P q X α ··· α n f l α ( τ ) ··· f l n α n ( τ n ) θ τ ··· τ n δδε ¯ l α ( τ ) ··· δδε ¯ l n α n ( τ n ) log D T e R t dτ P β [ ε + β ( τ ) φ − β ( τ )+ ε − β ( τ ) φ + β ( τ )] E (A6)where P q denotes all the permutation of the indexes k q ∈ { + , −} , such that there is a q number of minus super-operators. Performing the limit of ε ± α ( τ ) → one eventually obtains T log D T e − i R t dτ P α ( f + α ( τ ) φ − α ( τ )+ f − α ( τ ) φ + α ( τ )) E = ∞ X i =0 ( − i ) n Z t d τ n n X j =0 X α i , P q f l α ( τ ) ··· f l n α n ( τ n ) e C ¯ l ··· ¯ l n α ··· α n ( τ n ) , (A7)where the contribution given by the auxiliary degrees of freedom is described by the ordered cumulants: e C ¯ l ··· ¯ l n α ··· α n ( τ n ) = θ τ ··· τ n δ ε α ,τ ··· δ ε αn,τn log D T e R dτε α,τ φ α,τ E (cid:12)(cid:12)(cid:12) ε αi,τi =0 . (A8)A more explicit expression for the cumulants can be obtained by exploiting Ursell formula [31], that allows to write e C + ··· ¯ l n α , ··· α n ( τ n ) = X P ( | P |− − | − P | Y P ∈ P * T (cid:18) Y p ∈ P φ ¯ l p α p ( τ p )2 (cid:19)+ θ τ ...τ n (A9)where P is a partition of { , . . . , n } with cardinality | P | , P ∈ P is a set of the partition P , and p one element of theset P . Exploiting Eq. (A7) in Eq. (A1) and introducing the new symbol ˜ k n ( τ n ) = X α i , P q f l α ( τ ) ··· f l n α n ( τ n ) e C ¯ l ··· ¯ l n α ··· α n ( τ n ) , (A10)we eventually obtain Eq. (9) of the main text.We stress that the derivation provided in this Appendix holds also if φ ± in Eq. (A1) are stochastic processes and h· · ·i denotes the stochastic average. The final result is Eq. (A10) where e k → k and e C → C , i.e. Eq. (13). It isimportant to remark that in general e C and C differ because quantum trace and stochastic average have differentproperties. Appendix B
In this Appendix we provide calculations leading to Eq. (21). The map N t in Eq. (19) is TP if Tr[ N t ˆ ρ S ] = Tr[ˆ ρ S ] or equivalently ∂ t Tr[ N t ˆ ρ S ] = 0 . Performing the trace and taking the time derivative we get i∂ t Tr {N t [ˆ ρ S ] } =Tr ( T (cid:20) Z t d τ n ∞ X n =1 n X q =1 X α i , P q ( − i ) n f + α ( t ) f l α ( τ ) ··· f l n α n ( τ n ) C − ¯ l ··· ¯ l n α ··· α n ( τ n ) − iA + ( τ n ) ! × exp " ∞ X n =1 Z t d τ n (cid:16) ( − i ) n k n ( τ n ) − O n ( τ n ) (cid:17) ˆ ρ S ) . (B1)Exploiting the cyclicity of the trace we can rewrite the equation as i∂ t Tr {N t [ˆ ρ S ] } =Tr ( T (cid:20) Z t d τ n ∞ X n =1 n X q =1 X α i , P q ( − i ) n ˆ f α ( t ) f ↼ l α ( τ ) ··· f ↼ l n α n ( τ n ) C − ¯ l ··· ¯ l n α ··· α n ( τ n ) − i ˆ A ( τ n ) ! × exp " ∞ X n =1 Z t d τ n (cid:16) ( − i ) n k ↼ n ( τ n ) − O ↼ n ( τ n ) (cid:17) ˆ ρ S ) = 0; (B2)where the super-operator k ↼ ( τ n ) is k ( τ n ) where all the f l i have been replaced by f ↼ l i . From the above equation weimmediately obtain Eq. (21) for the operator ˆ A ( τ n ) . This condition is rather cumbersome, and one might argue thata simpler one can be obtained by expanding the operators ˆ O n ( τ n ) over the set { ˆ f α } as follows: ˆ O n ( τ n ) = X α i ˆ f α ( τ ) ··· ˆ f α n ( τ n ) K α ··· α n ( τ n ) , (B3)where the K α ··· α n are complex kernels. However, we now show that a decomposition of this kind works only in theGaussian case, while it is not instructive for noises with cumulants higher than the second. With the choice (B3), thestochastic map (17) generates the average dynamics (19) with (see Appendix C) O n ( τ n ) = n X q =0 X α i , P q f l α ( τ ) ··· f l n α n ( τ n ) K ⋄ α ··· α n ( τ n ) , (B4)where ⋄ = + when q is even, ⋄ = − when q is odd, and we have introduced K ± = ( K ± K ∗ ) / . One notices thatthe structure of the super-operators O n ( τ n ) is close to that of k n ( τ n ) . The relevant difference is that, unlike C , K depends only on the number q of plus signs of the array ( l , · · · l n ) , not on the array itself (note different superscripts).Substituting Eqs. (13) and (B4) in Eq. (19) one finds N t [ˆ ρ S ] = T exp " Z t d τ n ∞ X n =1 n X q =1 X α i , P q ( − i ) n f l α ( τ ) ··· f l n α n ( τ n ) (cid:16) C ¯ l ··· ¯ l n α ··· α n ( τ n ) − K ⋄ α ··· α n ( τ n ) (cid:17) ˆ ρ S . (B5)By performing the trace of this equation, one finds that the averaged map N t is TP only if the following conditionsare satisfied: ∀ n ∈ N , q ≤ n C − ¯ l ··· ¯ l n α ··· α n ( τ n ) − K + α ··· α n ( τ n ) = 0 , q even C − ¯ l ··· ¯ l n α ··· α n ( τ n ) − K − α ··· α n ( τ n ) = 0 , q odd . (B6)We recall that ¯ l = − because of the trace operation (see Eq. (14)). Since the identities of Eq. (B6) must hold forany array ( l , · · · l n ) , they represent a system of n − independent equations. However, at any order n there are onlytwo free variables: K + α ··· α n ( τ n ) and K − α ··· α n ( τ n ) . This is a consequence of the fact that C depends on the array ( − , ¯ l , · · · ¯ l n ) , while K depends only on the number of minus signs in such an array. Accordingly, we can concludethat the system (B6) admits solution only for n ≤ , and that the decomposition (B3) allows to obtain the conditionfor a TP map only in the Gaussian case. Appendix C
In this Appendix we derive Eqs. (B4) for the modified map N t . It is convenient to introduce the left-right formalismdenoting by a superscript L (R) the operators acting on ˆ ρ from the left (right) [29]. The map N t = h Ξ t ρ Ξ † t i , definedin the main text, can be then written as: N t = * T exp (cid:26) ∞ X n =1 ( − i ) n Z t d τ n O Ln ( τ n ) + ( − ) n O Rn ( τ n ) † − i Z t dτ X α [ f Lα ( τ ) φ Lα ( τ ) + f Rα ( τ ) φ Rα ( τ )] (cid:27)+ . (C1)We introduce the super-operator O n ( τ n ) : O n ( τ n ) = O Ln ( τ n ) + ( − ) n O Rn ( τ n ) † , (C2)which, making use of Eq. (B3) and of the identity K ± = ( K ± K ∗ ) / , can be rewritten as follows: O n ( τ n ) =[ f Lα ( τ ) ··· f Lα n ( τ n ) + ( − ) n f Rα ( τ ) ··· f Rα n ( τ n )] K + α , ··· ,α n ( τ n )+[ f Lα ( τ ) ··· f Lα n ( τ n ) − ( − ) n f Rα ( τ ) ··· f Rα n ( τ n )] K − α , ··· ,α n ( τ n ) . (C3)Further rewriting f R/L in terms of f ± = ( f L ± f R ) / , we obtain O n ( τ n ) = (cid:8) [ f + α ( τ ) + f − α ( τ )] ··· [ f + α n ( τ n ) + f − α n ( τ n )] + ( − ) n [ f + α ( τ ) − f − α ( τ )] ··· [ f + α n ( τ n ) − f − α n ( τ n ] (cid:9) K + α , ··· ,α n ( τ n )+ (cid:8) [ f + α ( τ ) + f − α ( τ )] ··· [ f + α n ( τ n ) + f − α n ( τ n )] − ( − ) n [ f + α ( τ ) − f − α ( τ )] ··· [ f + α n ( τ n ) − f − α n ( τ n ] (cid:9) K − α , ··· ,α n ( τ n )= n X j =1 X α i , P q (cid:26)(cid:2) f l α ( τ ) ··· f l n α n ( τ n )+( − ) q f l α ( τ ) ··· f l n α n ( τ n ) (cid:3) K + α , ··· ,α n ( τ n )+ (cid:2) f l α ( τ ) ··· f l n α n ( τ n ) − ( − ) q f l α ( τ ) ··· f l n α n ( τ n ) (cid:3) K − α , ··· ,α n ( τ n ) (cid:27) , (C4)where P q is the permutation of the indexes l i such that q is the number of plus signs in the array ( l , · · · , l n ) , and α i denotes the sum over all { α , · · · , α n } . The last line of this equation shows that the terms proportional to K + ( K − ) with odd (even) q are zero, allowing to rewrite the O n ( τ n ) as follows O n ( τ n ) = n X q =0 X α i , P q f l α ( τ ) · · · f l n α n ( τ n ) K ⋄ α ··· α n ( τ n ) , (C5)where ⋄ = + when q is even, ⋄ = − when q is odd, providing us with Eq. (B4). Appendix D
In this Appendix we provide explicit calculation for the series expansion of Eqs.(35)-(37). Since our aim is to providea series in powers of interaction operators ˆ f α , we need to rearrange the series in Eq. (33) because this contains terms0with powers of interaction operators in the range [ n, n ] . The strategy is to disentangle the two series in Eq. (33) byexploiting the properties of the Cauchy product of two series: ˆΞ t = T ∞ X n =0 ∞ X k =0 n ! k ! (cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n (cid:18) − Z t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) (cid:19) k . (D1)Notice that, at each increment of k the order of the ˆ f α increases by two, while at each increment of n it increases justby one. We replace k → ⌊ k/ ⌋ in order to have the same order of operators ˆ f α in both series and obtain: ˆΞ t = ∞ X n =0 ∞ X k =0 n ! ⌊ k/ ⌋ ! T (cid:26) (cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n (cid:18) − Z t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) (cid:19) ⌊ k/ ⌋ (cid:27) . (D2)The extra factor one half appearing in the second line of this equation is needed to compensate the extra contributionsgiven by the introduction of ⌊ k/ ⌋ , that counts twice the terms of the series ( ⌊ (2 k + 1) / ⌋ = ⌊ (2 k ) / ⌋ ). We can nowexploit again the properties of the Cauchy product in order to entangle back the series, obtaining ˆΞ t = ∞ X n =0 ( − i ) n ˆ ξ t (D3)where ˆ ξ t = n X k =0 n − ⌊ k/ ⌋ )! ⌊ k/ ⌋ ! T (cid:26) (cid:18)Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n − k (cid:18)Z t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) (cid:19) ⌊ k/ ⌋ (cid:27) . (D4)This series is such that the order n coincides with the power of operators ˆ f α . In order to make this fact more evidentwe replace k → k , to get ˆ ξ t = T ∞ X n =0 ⌊ n/ ⌋ X k =0 n − k )! k ! T (cid:26)(cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n − k (cid:18) − Z t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) (cid:19) k (cid:27) . (D5)This equation shows that at each increment of k one removes two terms of the type R t dτ ˆ f α ( τ ) φ α ( τ ) and adds oneterm of the type R t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) . In order to obtain Eq. (35) one only needs to make the time orderingexplicit by conditioning the integrals with unit step functions. Doing so we eventually obtain Eq. (35) in the maintext, i.e. ˆ ξ n = Z t d τ n X α i ˆ f α ( τ ) ··· ˆ f α n ( τ n ) ⌊ n/ ⌋ X k =1 X P ai K α α ( τ a , τ a ) ··· K α k − α k ( τ a k − , τ a k ) φ α k +1 ( τ a k +1 ) ··· φ α n ( τ a n ) θ τ a k +1 ··· τ an . (D6)In order to find Eq. (36), we replace Eq. (D1) in Eq. (32) obtaining h i ∂ t ˆΞ t − ˆ f α ( t ) φ α ( t ) ˆΞ t i = T (cid:26) − i Z t dτ ˆ f α ( t ) ˆ f β ( τ ) K αβ ( t, τ ) ∞ X n =0 ∞ X k =0 n ! k ! (cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n (cid:18) − Z t d τ ˆ f α ( τ ) ˆ f α ( τ ) K α α ( τ ) (cid:19) k (cid:27) . (D7)Following the steps that lead to Eq. (D3) one obtains Eq. (37). Appendix E
In this Appendix we show how one can obtain Eq. (44) by applying Wick’s theorem to Eq. (42). As preliminarysteps we separately study the following two cases: T ( ˆ f ( τ ) e − i R t dτ ˆ f ( τ ) φ ( τ ) ) = T ( ˆ f ( τ ) ∞ X n =0 n ! (cid:18) − i Z t dτ ˆ f ( τ ) φ ( τ ) (cid:19) n ) (E1)1and Z t dτ T n K ( t, τ ) ˆ f ( τ ) e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) o = T (Z t dτ K ( t, τ ) ˆ f ( τ ) ∞ X n =0 n ! (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n ) (E2)where for the sake of simplicity we have dropped the indexes α i . Let us start with Eq. (E1). Each term of the equationis composed by a string of n integrated operators: Z t d τ n T [ ˆ f ( τ ) ˆ f ( τ ) ··· ˆ f ( τ n )] . (E3)Applying Wick’s theorem one obtains: Z t d τ n T [ ˆ f ( τ ) ˆ f ( τ ) ··· ˆ f ( τ n )] = Z t d τ n ˆ f ( τ ) T [ ˆ f ( τ ) ··· ˆ f ( τ n )] + n X i =1 Z t d τ n C ( τ , τ i ) T [ ˆ f ( τ ) ··· ˆ f ( τ i − ) , ˆ f ( τ i +1 ) ··· ˆ f ( τ n )]= Z t d τ n ˆ f ( τ ) T [ ˆ f ( τ ) ··· ˆ f ( τ n )] + n Z t d τ n C ( τ , τ ) T [ ˆ f ( τ ) ··· ˆ f ( τ n )] , (E4)where C ( τ , τ ) ≡ [ ˆ f ( τ ) , ˆ f ( τ )] θ τ ,τ is a Wick contraction. Exploiting Eq. (E4) in Eq. (E1) we are able to extract theoperator ˆ f ( τ ) from the time ordering, at the cost of introducing an extra term proportional to the contraction: T ( ˆ f ( τ ) e − i R t dτ ˆ f α ( τ ) φ α ( τ ) ) = (cid:18) ˆ f ( τ ) − i Z t dτ C ( τ , τ ) (cid:19) ∞ X n =0 n ! T (cid:18) − i Z t dτ ˆ f α ( τ ) φ α ( τ ) (cid:19) n = (cid:18) ˆ f ( τ ) − i Z t dτ C ( τ , τ ) (cid:19) T e − i R t dτ ˆ f α ( τ ) φ α ( τ ) (E5)The next step is to show how one can extract the term ˆ f ( τ ) from the time ordering in Eq. (E2). Applying Wick’stheorem for the n -th term of the series one gets: T (Z t dτ K ( t, τ ) ˆ f ( τ ) 1 n ! (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n ) = Z t dτ K ( t, τ ) ˆ f ( τ ) T ( n ! (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n ) − Z t dτ K ( t, τ ) T (Z t d τ [ C ( τ , τ ) ˆ f ( τ ) + C ( τ , τ ) ˆ f ( τ )] K ( τ ) 1( n − (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n − ) = Z t dτ K ( t, τ ) ˆ f ( τ ) T ( n ! (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n ) − Z t dτ Z t d τ K ( t, τ ) C ( τ , τ ) T ( K ( τ , τ ) ˆ f ( τ ) 1( n − (cid:18) − Z t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:19) n − ) (E6)where K ( τ , τ ) = K ( τ , τ ) + K ( τ , τ ) . Applying this result to Eq. (E2) one obtains: T (cid:26) ˆ f ( t ) Z t dτ K ( t, τ ) ˆ f ( τ ) e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:27) = −T (cid:26)Z t dτ ˆ f ( τ ) K (2) ( t, τ ) e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:27) + ˆ f ( t ) Z t dτ K ( t, τ ) (cid:18) ˆ f ( τ ) − i Z t dτ C ( τ , τ ) φ ( τ ) (cid:19) T e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) . (E7)As one can see, unlike the previous case, in this specific case we are not able to close the equation at the first stepbecause we obtain an extra term. However, such a term has the same structure as the one in the l.h.s. of Eq. (E7),allowing to re-apply this equation to it. Iterating this procedure we obtain T n K ( t, τ ) ˆ f ( τ ) e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) o = ∞ X n =1 Z t dτ K ( n ) ( t, τ ) ˆ f ( τ ) T e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) − T (cid:26)Z t dτ ˆ f ( τ ) K ( ∞ ) ( t, τ ) e − R t d τ ˆ f ( τ ) ˆ f ( τ ) K ( τ ) (cid:27) . (E8)2with K (1) = K and K ( n ) ( t, τ ) = Z t d τ K ( n − ( t, τ ) C ( τ , τ ) K ( τ , τ ) . (E9)Eq. (E8) shows that one can extract the term ˆ f ( τ ) from the time ordering only if lim n →∞ K ( n ) = 0 (E10)(otherwise an extra term would be left, not allowing to close the procedure). This condition can be understood as aconstraint on the convergence of the perturbative series. With Eqs. (E5) and (E8) in hand it is straightforward toobtain Eq. (44) from Eq. (43). [1] N. Gisin, Phys. Rev. Lett. ,1657 (1984); L. Diósi, Phys.Lett. A , 451 (1986), .[2] A. Bassi, G.C. Ghirardi, Phys. Rep. , 257 (2003); A.Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht,Rev. Mod. Phys. , 471 (2013).[3] V. P. Belavkin, Phys. Lett. A , 355 (1989); J. Math.Phys. , 2930 (1990).[4] N. G. van Kampen, Stochastic Processes in Physics andChemistry , (North Holland, 2006).[5] A. Barchielli, M. Gregoratti.
Quantum trajectories andmeasurements in continuous time: the diffusive case ,Lect. Notes Phys. Vol. 782. (Springer, 2009).[6] A. Smirne, H.-P Breuer, J. Piilo, B. Vacchini, Phys. Rev.A , 062114 (2010); T. Guerin, O. Benichou, and R.Voituriez, Nat. Chem. , 568 (2012); S. Gröblacher, A.Trubarov, N. Prigge, G.D. Cole, M. Aspelmeyer, and J.Eisert, Nat. Commun. , 7606 (2015).[7] I. de Vega, D. Alonso, Rev. Mod. Phys. , 015001(2017).[8] L. Diósi, W. T. Strunz, Phys. Lett. A , 569 (1997).[9] S. L. Adler, A. Bassi, J. Phys. A ibid. , 012116 (2009).[11] L. Ferialdi and A. Bassi Phys. Rev. A , 022108 (2012).[12] A. Smirne, B. Vacchini, A. Bassi, Phys. Rev. A ,062135 (2014); A. Smirne, A. Bassi, Sci. Rep. , 12518(2015); M. Toros, G. Gasbarri, A. Bassi, Phys. Lett. A , 3921 (2017).[13] J. Gambetta, H. M. Wiseman, Phys. Rev. A ,012108(2002); L. Diósi, Phys. Rev. Lett. , 080401 (2008);H.M. Wiseman and J.M. Gambetta, Phys. Rev. Lett. , 140401 (2008); L. Diósi, Phys. Rev. Lett. ,149902(E) (2008).[14] A. Barchielli, V.P. Belavkin, J. Phys. A: Math. Gen.
293 (1995); A. Barchielli, A.M. Paganoni, F.Zucca, Stoch. Proc. Appl.
69 (1998); A. Barchielli,M. Gregoratti, Phil. Trans. R. Soc. A , 5364 (2012);A. Barchielli, C. Pellegrini, F. Petruccione, Phys. Rev. A , 063814 (2012); I. Semina, V. Semin, F. Petruccione,A. Barchielli, Open Syst. Inf. Dyn. , 1440008 (2014).[15] W. T. Strunz and T. Yu, Phys. Rev. A , 52115 (2004);J. Jing and T. Yu, Phys. Rev. Lett. , 240403 (2010);X. Zhao, J. Jing, B. Corn, and T. Yu, Phys. Rev. A ,32101(2011).[16] L. Ferialdi, A. Bassi, Europhys. Lett. , 30009 (2012).[17] A. Tilloy, Quantum , 29 (2017).[18] L. Diósi, L. Ferialdi, Phys. Rev. Lett. , 200403 (2014).[19] A. A. Budini, Phys. Rev. A , 052101 (2015).[20] L. Diósi, Phys. Rev. A , 034101 (2012); N. Megier,W. T. Strunz, C. Viviescas, K. Luoma, arXiv:1710.08359(2017).[21] M.G. Genoni and M. G. A. Paris, Phys. Rev. A ,052341 (2010); C. Benedetti, F. Buscemi, P. Bordone,and M. G. A. Paris, Phys. Rev. A , 032114 (2014); N.Leigh, G. A. Paz-Silva, L. Viola, Phys. Rev. Lett. ,150503 (2016).[22] W. F. Stinespring, Proc. Am. Math. Soc , 211 (1955).[23] A. Bassi, E. Ippoliti, B. Vacchini, J. Phys. A , 8017(2005).[24] E. Engel, R. M. Dreizler, Density Functional Theory (Springer-Verlag Berlin Heidelberg, 2011).[25] K. Furutsu, Nat. Bur. (1964); E.A. Novikov, J. Eksp.Teor. Phys. , 5 (1964); K. Sobczyk, Stochastic WavePropagation (Elsevier, Amsterdam 1985).[26] A. Budini, Phys. Rev. A , 012106 (2000).[27] G. Gasbarri, L. Ferialdi, Phys. Rev. A , 022114 (2018).[28] L. Ferialdi, Phys. Rev. Lett. , 120402 (2016); Phys.Rev. A , 052109 (2017).[29] K. Chou, Z. Su, B. Hao, and L. Yu, Phys. Rep. , 1(1985); L. Diosi, Found. Phys. , 63 (1990); L. Diosi,Physica A , 517 (1993).[30] M. Reed, B. Simon, Functional Analysis , (AcademicPress, 1981).[31] H. D. Ursell, Math. Proc. Cam. Phil. Soc. , CambridgeUniversity Press, (1927); S. Garrett, Comm. Math. Phys.42