A discrete memory-kernel for multi-time correlations in non-Markovian quantum processes
AA discrete memory-kernel for multi-time correlations in non-Markovian quantum processes
Mathias R. Jørgensen ∗ and Felix A. Pollock † Department of Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark School of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia (Dated: July 8, 2020)Efficient simulations of the dynamics of open systems is of wide importance for quantum science and tech-nology. Here, we introduce a generalization of the transfer-tensor , or discrete-time memory kernel, formalismto multi-time measurement scenarios. The transfer-tensor method sets out to compute the state of an open few-body quantum system at long times, given that only short-time system trajectories are available. Here, we showthat the transfer-tensor method can be extended to processes which include multiple interrogations (e.g. mea-surements) of the open system dynamics as it evolves, allowing us to propagate high order short-time correlationfunctions to later times, without further recourse to the underlying system-environment evolution. Our approachexploits the process-tensor description of open quantum processes to represent and propagate the dynamics interms of an object from which any multi-time correlation can be extracted. As an illustration of the utility of themethod, we study the build-up of system-environment correlations in the paradigmatic spin-boson model, andcompute steady-state emission spectra, taking fully into account system-environment correlations present in thesteady state.
I. INTRODUCTION
Modelling the dynamical properties of open quantum sys-tems is an outstanding challenge in modern physics. Ap-plications range from studying impurity dynamics in ultra-cold quantum gases or strongly correlated materials [1–3],to modelling emission properties of organic molecules or ar-tificial atoms [4, 5]. While a plethora of techniques existfor the simulation of the time-dependence of observables ata single point in time [6, 7], usually encoded in the time-dependent state of an open system, there are relatively fewthat are aimed at computing multi-time correlation functionsin full generality. These correlation functions are critical tomodelling spectral responses, such as those measured in ul-trafast spectroscopy experiments [8, 9], used to directly in-fer quantum dynamical behaviour in, for example, photosyn-thetic processes [10].If the coupling between the system under study and its sur-rounding environment is sufficiently weak, and the reorga-nization timescale of the environment sufficiently short, thatmemory effects can be neglected, then multi-time correlationscan be expressed in terms of two-time correlations throughthe so-called quantum regression theorem [11]. Under an ad-ditional stationarity assumption, these two-time correlationsare themselves expressible in terms of single-time observableevolution [6]. However, for more general non-Markovian dy-namics, environment-mediated memory effects preclude suchsimplifications, and correlations across multiple times are notfully determined by those across fewer [12–14].Accurately computing these correlations typically requiressophisticated numerical methods, such as those based on pathintegrals [15, 16] and quantum Monte Carlo [17, 18], andother methods that rely on simplifying details of the modelin question [19]. In many cases, these methods suffer frompoor convergence with system size, with system-environment ∗ [email protected] † [email protected] interaction strength, and with the total simulation time. Thislimits the ability to apply them to important processes involv-ing large, complex systems or where long-time dynamics isimportant.For (single-time) state evolution, these kinds of scaling dif-ficulties can, in some cases, be mitigated with the aid of a memory kernel [6, 20, 21], a mapping that quantifies howa system’s past trajectory contributes to its future evolution.Once this is determined, the system’s state can be propagatedin time with an efficiency that does not depend on the detailsof the underlying system-environment dynamics [22]. Thisis perhaps most clearly evident in the discrete-time mem-ory kernel approach, known as the transfer-tensor method(TTM) [23, 24] (closely related to an earlier approach spe-cific to the spin-boson model [25]). The TTM can be treatedas a black-box, taking short-time non-Markovian dynamicalmaps, computed using some other technique, as input, andself-consistently propagating them to arbitrarily long times.It has been used to great success in simulating properties ofsystems in contact with a bosonic bath [26, 27]; however, it isnot directly applicable to the computation of multi-time cor-relations.In this paper, we develop a generalization of the transfer-tensor formalism to a scenario involving sequential measure-ments of an non-Markovian open system as it evolves, result-ing in an efficient discrete-time memory kernel method for thepropagation of multi-time correlation functions. In contrastto previous multi-time memory kernel methods, for whichcorrelations between each set of operators require a sepa-rate kernel [28], our approach utilizes the powerful process-tensor formalism to encode all correlations of a given orderinto a single positive operator [16, 29]. As an illustration ofthe utility of our generalized method, we study the build-upof system-environment correlations in the paradigmatic spin-boson model, and compute steady-state emission spectra, tak-ing into account system-environment correlations present inthe steady state. a r X i v : . [ qu a n t - ph ] J u l II. BACKGROUND
In this section we introduce the basic notions of open quan-tum dynamics and multi-time correlation functions. Often,one is interested in the probability distribution describing ob-served outcomes in a sequential measurement scenario, andwe discuss how such a probability distribution can be conve-niently expressed in terms of a process tensor [29]. Lastlywe outline the standard transfer-tensor method, as it appliesto the propagation of the system’s state.
A. Open dynamics and correlation functions
We are here concerned with the simulation of a stochasticprocess undergone by an open quantum system over the setof discretized times { t n = nδt : n ∈ N } , where δt definesthe temporal resolution of the simulation . In a general quan-tum stochastic process of the sort we consider here, an opensystem ( S ) evolves along with its environment ( E ), while si-multaneously being monitored and manipulated by a putativeexperimenter. Its dynamics therefore has two contributions:the first is the dynamical evolution the open system wouldundergo between consecutive discrete times if it were not ex-ternally perturbed, and the second is the influence of mea-surements and other external interventions. Here, we imag-ine the latter are implemented sufficiently fast that they can beconsidered instantaneous with respect to δt ; continuous ma-nipulations can be approximated in a controlled fashion byreducing δt appropriately [16, 30].The free dynamical evolution of an open system can alwaysbe modelled by considering that the joint SE dynamics canbe described by a time evolution superoperator Λ t,s governedby the Liouville-von Neumann equation, ∂ t Λ t,s = L t Λ t,s . (1)The time evolution is generated by the Liouvillian superoper-ator defined by L t ( • ) ≡ − i [ H ( t ) , • ] , where the Hamiltonian H describes the system, the environment and their interac-tion. With the initial condition Λ s,s = I (the identity super-operator), this has the formal solution Λ t,s = T ← e (cid:82) ts dr L r ,with T ← indicating time ordering of the following exponen-tial. We will henceforth specialise to the case where the SE Liouvillian is time-independent: L t = L ∀ t , such that, withrespect to our discrete time grid, we can write Λ n := Λ t n , = e nδt L . In this case, we have that Π n = e nδt L (Π ) , (2)where Π n denotes the full system-environment state at time t n , and the above equation relates this state to the full ini-tial state Π . Throughout we focus on the scenario in which For the purposes of the transfer-tensor method we derive in later sections,this resolution need not have any relation to the intrinsic timescales of theprocess to be simulated. However, many techniques for computing theinitial dynamics rely on δt being sufficiently small [7]. the dynamics is generated by a time-independent Hamilto-nian, however the results we present can straightforwardlybe generalized to a periodically time-dependent Hamiltonianfollowing the approach of Pollock et al. [24]. Furthermore wework in units in which (cid:126) = 1 , such that time is measured inunits of inverse energy.It is often the case, that we want to compute opera-tor expectation values or multi-time correlation functions.Suppose, that the system has an observable represented bythe Hermitian operator A with spectral resolution A = (cid:80) k a k | a k (cid:105)(cid:104) a k | , then the expectation value of the operatorat time t n is given in the Heisenberg picture by (cid:104) A ( t n ) (cid:105) = tr SE (cid:2) e it n H Ae − it n H Π (cid:3) = (cid:88) k a k tr SE (cid:2) | a k (cid:105)(cid:104) a k | e t n L (Π ) (cid:3) = (cid:88) k a k P ( a k , t n ) , (3)where we have introduced P ( a k , t n ) as the probability ofmeasuring the value a k of the observable A at time t n . Anoperator expectation value can be seen as a process in whichthe system evolves up to a time t n , and where an observable A is then projectively measured. Suppose now the system has asecond observable B = (cid:80) j b j | b j (cid:105)(cid:104) b j | , we can then considertwo-time correlation functions of the form (cid:104) A ( t n ) B ( t m ) (cid:105) = (cid:88) kj a k b j W ( a k , t n ; b j , t m ) , (4)where in this case we have introduced the function W (noticethat this is not a probability distribution) by W ( a k , t n ; b j , t m ) =tr SE (cid:104) | a k (cid:105)(cid:104) a k | e ( t n − t m ) L ( | b j (cid:105)(cid:104) b j | e t m L (Π )) (cid:105) . (5)Two-time correlation functions of this type do not corresponddirectly to any realizable sequence of measurements, and it isnot generally possible to establish a direct relation to an un-derlying unique probability distribution (In general, you canalways write these in terms of a linear combination of proba-bilities, albeit not from measurements in a single basis [31]).In spite of this, correlation functions find wide application.For our purposes the crucial thing is that in both cases thefunction we want to compute can be expressed as unitary evo-lution on the joint system-environment space, punctuated byan operation on the system.The formalism developed below focus on the case wherethe implemented operations corresponds to generalized mea-surements, this approach has the advantage that one can workstrictly with positive operators. The developed method, how-ever, can equally well be applied to more general correlationfunctions, such as the one considered above. B. General measurements and the process tensor
Generally, an instantaneous measurement on the open sys-tem is represented mathematically by a completely-positivesuperoperator O Sx n , which operates only on the system andnot on the environment. Here the index x n labels the set ofmeasurement outcomes of the specific measurement imple-mented at time t n . For a measurement implemented instanta-neously the full state transforms as Π n ; x n = O Sx n ⊗ I E (Π n ) . (6)In the remainder of the paper we will for convenience notwrite the identity map on the environment or the superscript S on the system superoperator explicitly. Performing a mea-surement is a non-deterministic operation, this means that theresulting state is generally conditioned on the measurementoutcome (indicated by the subscript), and sub-normalized(we take this to be implied by the presence of the sub-script). Throughout we will keep to the operational languagein which the considered operations correspond to a gener-alized measurement. However, the results we present canstraightforwardly accommodate superoperators representingmore abstract transformations useful in the computation ofe.g. nonequilibrium Green functions.Now, consider the statistics obtained by performing a mea-surement at every discrete timestep. The probability of realiz-ing a particular measurement record x n :0 = { x n , ..., x } thentakes the form P ( x n :0 ) = tr (cid:2) O x n e δt L O x n − ... e δt L O x (Π ) (cid:3) . (7)Notice that even though we are including a measurement ateach timestep, we are not in the final analysis forced to actu-ally implement a measurement, since the identity map corre-sponds to a trivial measurementLater when we present our generalized transfer-tensor for-malism, we will find it useful to note that the above proba-bility can be rewritten by means of the Choi-Jamiolkowskiisomorphism [29]. Suppose that we have available a collec-tion of copies of our system. For reasons which will be clearshortly we semantically divide the set of copies into inputsystems R j and output systems S j . Given this collection ofancilla systems we introduce the process tensor in its many-body Choi state representation [29] Υ n :0 = tr E (cid:2) e δt L n ... e δt L (cid:0) Ψ S n R n − ⊗ ... Ψ S R ⊗ Π S E (cid:1)(cid:3) . (8)Here the Liouville operator L j operates on the output space S j and the environment E , and furthermore Π S E denotes thearbitrary initial state of the environment and output space S .In the above Ψ SR ≡ (cid:80) da,b =1 | aa (cid:105)(cid:104) bb | denotes the unnormal-ized maximally entangled state of ancillary systems S and R ,where d is the Hilbert space dimension. Although the aboveconstitutes the Choi state representation of the process tensor,we will for simplicity refer to it as the process tensor through-out. Furthermore we define the Choi state of the measurementsequence by O x n :0 = O x n ⊗ I ⊗ O x n − ... ⊗ I ⊗ O x (cid:0) S n ⊗ Ψ R n − S n − ⊗ ... Ψ R S (cid:1) , (9)where O x j operates on output space S j . Throughout we willrefer to this simply as the measurement sequence. In terms of the process tensor and the measurement sequence, we can re-cast the probability of realizing the measurement record x n :0 in the form of a spatio-temporal Born rule [29] P ( x n :0 ) = tr (cid:2) O Tx n :0 Υ n :0 (cid:3) , (10)where the trace is over all input and output system spaces.Given this expression for the probability, it is clear that theprocess tensor constitutes a generalization of the open sys-tem state to multi-time measurement scenarios. Crucially forour purposes here, the process tensor makes it possible to de-velop a generalization of the transfer-tensor formalism, with-out having to deal with the complications of working withstates conditioned on particular measurement outcomes. C. Transfer-tensor formalism for quantum state evolution
Lastly we introduce the standard transfer-tensor formalism.Consider the simplest process in which temporal correlationsplay a role, namely that of a two-time measurement. Sucha process consists of an initial measurement (we call thisthe preparation procedure) at time t , followed by a secondmeasurement at a later time t n . In this case, the probabilityof realizing the process defined by the measurement record { x n , x } takes the form P ( x n , x ) = tr (cid:2) O x n e nδt L O x (Π ) (cid:3) . (11)where the maps are applied consecutively. To compute thisprobability, we must generally solve the following simulationproblem: First the initial system-environment state must besubjected to the preparation procedure, and then the resultingconditional state must be evolved up to timestep n accordingto the Liouville-von Neumann equation. The statistics associ-ated with the second measurement at timestep n , is then fullycharacterized by the conditional open system state ρ n ; x = tr E (cid:2) e nδt L O x (Π ) (cid:3) . (12)Computing the open system state by simulating the evolutionof the full system-environment state, generally constitutes ahighly inefficient computational task. For this reason it isdesirable to develop a dynamical description of the processexclusively within the open system subspace, which accountsfor the effects of the environment on the system dynamics,but otherwise eliminates all information on the environment.A first step towards such a reduced description, can betaken by pointing out that any given system-environment statecan be decomposed into an uncorrelated part and a correctionterm. We take the uncorrelated part to be given by the tensorproduct of the open system state and an environment refer-ence state τ E which we take to be independent of time. Thecorrection term accounts for any system-environment corre-lations, and for any discrepancy between the reduced envi-ronment state and the environment reference state. It thenfollows that we can write the initial system-environment stateas Π = ρ ⊗ τ E + χ , where χ is the correction term.Throughout we will refer to the uncorrelated first term as theproduct state projection. If we substitute this decompositioninto the right-hand side of Eq. (12), then we find that the opensystem state at time t n consists of a homogeneous and an in-homogeneous contribution: ρ n ; x = E n ( ρ x ) + J n ; x . (13)The inhomogeneous term J n ; x accounts for the discrepancybetween the actual initial system-environment state and theproduct state projection. We will discuss the problem of ac-counting for this term towards the end of this section. Thehomogeneous contribution is expressed in terms of the dy-namical map defined by E n ( • ) = tr E (cid:2) e nδt L ( • ⊗ τ E ) (cid:3) . (14)It relates the initial system state to the system state at alater time, under the assumption that the initial system-environment state is described by the product state projection.Among many possible definitions [32–35], the presence oftemporal correlations in the homogeneous open system dy-namics can be related to a violation of the semi-group prop-erty E n = E n . A violation of this property implies that thedynamical maps at different times must be computed individ-ually, which is computationally inefficient. Progress on thisproblem, can be made by noticing that the set of dynamicalmaps has an equivalent representation as a set of transfer ten-sors [23]. The transformation rule from one set to the other isthe following recursive relation: T n = E n − n − (cid:88) j =1 T j E n − j . (15)The homogeneous open system dynamics can equivalently beexpressed either in terms of the dynamical maps or the trans-fer tensors. In fact the transfer tensors can be seen as the dis-cretized analogue of the Nakajima-Zwanzig memory kernel[23, 24].The feature of the transfer-tensor representation that makesit interesting can be appreciated as follows: Notice that while T = E , the two-step transfer tensor represents the deviationbetween the two-step dynamical map and the dynamical mapcomposed of consecutive applications of the one-step map T = E − E E . Generally the transfer tensor T n encodesthe contribution to the dynamical map E n from temporal cor-relations in the dynamics extending for a time nδt . For anopen system dynamics which has a bounded memory time,we thus expect that the contribution from the transfer tensor T n to the open system dynamics should approach zero as n increases. If this is approximately the case, then the recursiverelation Eq. (15) can be used to iteratively construct approx-imate dynamical maps extending over arbitrarily long times,given only a finite set of transfer tensors. We postpone a dis-cussion of the approximation error to the next section, whereit can be approached for the multi-time case. Before pro-ceeding we point out that the transfer-tensor representation,does not free us from having to construct the set of dynam-ical maps up to a time comparable to the correlation time ofthe environment. Rather they provide an efficient represen-tation of the information contained in the constructed maps,and allows us to efficiently obtain dynamical maps extendingover longer times than initially simulated, if the open systemdynamics is finitely correlated in time. Lastly we must return to the inhomogeneous term. Theproblem of accounting for a non-zero inhomogeneous termhas been discussed by Buser et al. [27], who proposed the fol-lowing strategy: (i) The dynamical maps are generated froma simulation method of choice, and the transfer tensor are re-constructed. (ii) The reduced state corresponding to corre-lated initial conditions is simulated using a numerical methodcapable of dealing with correlated initial conditions. (iii) Theinhomogeneous term is inferred by looking at the differencebetween the simulated open system state and the homoge-neous contribution. If it is found that the inhomogeneousterm vanishes at sufficiently long times, then the long-timebehaviour can be studied exclusively in terms of the homoge-neous contribution. III. MULTI-TIME MEASUREMENTS
In the last section we discussed the standard transfer-tensormethod for quantum state evolution, and described how theprocess tensor is the natural generalization of the quantumstate to multi-time measurement scenarios. In this sectionwe make use of the process tensor, to construct a generalizedtransfer-tensor formalism, which is applicable to any multi-measurement process and not only to quantum state evolu-tion. In addition we analyse the error growth associated withimplementing a memory time cutoff, and provide an upperbound on the error.
A. Multi-time transfer tensors
Similarly to the two-time measurement case, we write theinitial system-environment state as the sum of a product stateprojection and a correction term. The process tensor then de-composes into a homogeneous and an inhomogeneous contri-bution: Υ n :0 = E n :0 ⊗ ρ + J n :0 . (16)The inhomogeneous term J n :0 accounts for any differencesbetween the actual initial system-environment state and theproduct-state projection. The homogeneous term is given bythe tensor product of the initial system state and the dynami-cal tensor defined by E n : k ≡ tr E (cid:2) e δt L n e δt L n − ... e δt L k +1 (cid:0) Ψ S n R n − ⊗ ... Ψ S k +1 R k ⊗ τ E (cid:1)(cid:3) . (17)The dynamical tensor provides the multi-time generalizationof the dynamical maps introduced for the two-time measure-ment (see Eq. (14)), and conceptually we can think of the dy-namical tensor as the Choi state representation of a sequenceof correlated dynamical maps on the open system space (seeFig. 1).We now derive a generalized transfer-tensor representa-tion of the dynamical tensor, by adapting the approach devel-oped for the two-time transfer-tensor formalism by Pollock etal. [24]. First we define the following projection superopera-tor acting on the dynamical tensor P j E n : k ≡ E n : j ⊗ E j : k for k < j < n (18)and its complement Q j ≡ I − P j . These projection opera-tors are analogous to the standard Nakajima-Zwanzig projec-tors [6], notice, however, that they are defined directly on thedynamical tensor and not on the underlying SE dynamics .Making use of the identity, we can iteratively decompose thedynamical tensor as E n : k = ( P n − + Q n − ( P n − + Q n − ( ... ))) E n : k = n − (cid:88) j = k +1 T n : j ⊗ E j : k + T n : k (19)where we have introduced the set of generalized transfer ten-sors T n : k . It is straightforward to show that these are relatedto the dynamical tensor by the recursive relation T n : k = Q n − ... Q k +1 E n : k = E n : k − n − (cid:88) j = k +1 T n : j ⊗ E j : k . (20)Notice that similarly to the standard case we have for the one-step transfer tensor T = E , while for the two-step trans-fer tensor T = E − E ⊗ E . Generally the trans-fer tensor T n :0 quantifies the contribution to the dynamicaltensor E n :0 from temporal correlations extending for a timeduration nδt . In contrast to the standard case, however, theformulation here allows for experimental interventions at allintermediate times. B. Long-time propagation and error growth
Now suppose that a simulation of the dynamical tensor E n :0 can be performed up to a timestep l . Then using the above re-cursive relation we can iteratively reconstruct the set of trans-fer tensors T l : k for all k < l . This is possible since the openevolution under study is generated by a time-independentHamiltonian, which implies that the dynamical tensors E n : k ,and consequently the transfer tensors, become translation in-variant in the sense that E n : k and E n − k :0 are equivalent up toa translation of the system input-output spaces on which theyare defined. For a dynamics finitely correlated in time, we ex-pect that the transfer tensors describing temporal correlationsextending for longer than the environment correlation time,to contribute negligibly to the dynamics of the system. If weassume that the timestep l is sufficiently large in this sense, Notice that due to the causal nature of the dynamical tensor, one can obtain E j : k from E n : k by taking the trace over all input and output spaces fromtimestep n down to the output space at timestep j . In the dilated/underlying picture, the action of the projection operators canbe equivalently represented by the usual Nakajima-Zwanzig projectors
Figure 1. Diagram illustrating the relation between the multi-timedynamical maps and the transfer tensors. The transfer tensors de-compose the dynamical maps into elements encoding memory ef-fects across a specific timescale. The generalized approach includestransfer tensors propagating memory effects across an implementedmeasurement. then to a good approximation the homogeneous process canbe modelled by the set of truncated dynamical maps E ( l ) n :0 = l (cid:88) k =1 T n : n − k ⊗ E ( l ) n − k :0 for n > l. (21)The crucial thing here is that given only the initially simulatedset of dynamical tensors, and the reconstructed set of trans-fer tensors, we can compute approximate dynamical tensorsextending over longer times than that initially simulated.To quantify the error induced by truncating the transfer ten-sor expansion, we use the Frobenius norm of the differencebetween the exact and the reconstructed dynamical tensor.Writing out the individual terms in the difference, and mak-ing use of the triangle inequality and the Cauchy-Schwartzinequality, we straightforwardly find that the error is upperbounded as || E n :0 − E ( l ) n :0 || F ≤ n (cid:88) k = l +1 ( n + 1 − k ) || T k :0 || F . (22)If we take it as a condition for the TTM to be applicablethat the transfer tensor norm is monotonically decreasing fortimes longer than the imposed memory cutoff, then the errorgrowth is in the worst case quadratic, with a proportionalityconstant given by the leading transfer tensor norm (cid:107)T l +1:0 (cid:107) F .In many cases this error growth is sufficiently slow to enableaccurate predictions of the system steady-state. C. Generating closed dynamical relations
We now show how the standard transfer-tensor formalismcan be recovered from the general formulae. First we con-sider a measurement sequence consisting of a measurementat time t and one again at time t n . All other measurementsare taken to be identity maps, that is O x j = I for all j except j = 0 and j = n . Then we can take the trace in Eq. (10)over all intermediate output-input spaces, except for outputspace S n and the input-output space R S . In carrying outthe trace we are projecting an output-input space R j S j of theprocess tensor onto a maximally entangled state Ψ R j S j . FromEq. (16) this results in the two-time process tensor Υ n, = E n, ⊗ ρ + J n, , (23)where the inhomogeneous terms is obtained from J n :0 byprojections at intermediate times. The homogeneous contri-bution is in this case given in terms of the two-time dynamicaltensor obtained by subjecting Eq. (17) to the projection pro-cedure E n,k = tr E (cid:104) e ( n − k ) δt L n (Ψ S n R k ⊗ τ E ) (cid:105) . (24)The equivalent two-time transfer-tensor representation is ob-tained by subjecting the general relation Eq. (20) to the sameprojection onto maximally entangled states of intermediaryoutput-input spaces, that is T n,m = E n,m − n − (cid:88) k = m +1 T n,k ∗ E k,m . (25)For convenience, we have introduced the star notation to de-note the projection of the intermediate boundary output-inputspace, that is T n,k ∗ E k,m ≡ tr R k S k [Ψ R k S k ( T n,k ⊗ E k,m )] . (26)The above results show, that in considering the two-time mea-surement, we recover the standard transfer-tensor formalism,with the trivial difference that instead of relations involvingcompositions of superoperators, we have a projection of theintermediate output-input space as a consequence of workingwithin the Choi state representation.In addition to recovering the standard two-time measure-ment results, the generalized formalism makes it possible tostudy three-time measurement statistics. We consider a mea-surement sequence consisting of an initial measurement attime t followed by a measurement at time t m and one againlater at time t n . All other measurements are taken to beidentity maps. Then we take the trace in Eq. (10) over alloutput-input spaces except for R S , R m S m and S n . Tak-ing the trace corresponds to projecting intermediate output-input spaces onto maximally entangled states, and similarlyto above this results in the three-time process tensor Υ n,m, = E n,m, ⊗ ρ + J n,m, , (27)where the inhomogeneous terms is obtained from J n :0 byprojections at all intermediate times except for time t m . The homogeneous term is expressed in terms of the three-time dy-namical tensor obtained by subjecting Eq. (17) to the projec-tion procedure E n,m,k = tr E [ e ( n − m ) δt L n e ( m − k ) δt L m (Ψ S n R m ⊗ Ψ S m R k ⊗ τ E )] . (28)Furthermore, by subjecting Eq. (20) to the projection proce-dure, it follows that the three-time dynamical tensor has thefollowing transfer-tensor representation T n,m,k = E n,m,k − m − (cid:88) j = k +1 T n,m,j ∗ E j,k − T n,m ⊗ E m,k − n − (cid:88) j = m +1 T n,j E j,m,k . (29)Notice that we obtain a closed dynamical relation, in the sensethat the three-time dynamical tensors can be decomposed intocontributions from the two-time dynamical tensors, the two-time transfer tensors and the three-time transfer tensors. Therelation between the transfer tensors and the dynamical tensoris illustrated in Fig. 1.Furthermore, the three-time transfer tensors have the fea-ture that if we want to compute steady-state two-point cor-relation functions, we can simply specify a sufficiently largevalue of m in the above equation. The equation is closed inthe sense that the three-time transfer tensors and the three-time dynamical maps, need not be computed for any othervalue of m . This makes it possible to compute steady-statecorrelation functions by first propagating the system into thesteady state, by a two-time transfer tensors propagation fora sufficiently long-time. Then an operation is performed onthe system (or an arbitrary operator is implemented). Fol-lowing the operation the system can be further propagatedusing the two- and three-time transfer tensors. Crucially,this can be done while taking full account of the system-environment correlations present in the steady state, some-thing which has not generally been possible before. Whensuch calculations are performed, the quantum regression the-orem is often applied [6], which corresponds to discarding thethree-time transfer-tensors.Lastly we mention the flexibility inherent in the transfer-tensor method in choosing a suitable temporal resolution.That is, even though we generate dynamical maps over atemporal grid with resolution δt , we can define the transfertensors with respect to a coarser grid, sδt for s some posi-tive integer, without increasing the overall simulation error.This means we can tailor the transfer-tensor decompositionto any desired temporal resolution of the system dynamics. Ifthe desired resolution has the feature that dynamics is well-captured by the single leading transfer tensor, then we recovera Markovian description of the process. IV. APPLICATION TO THE SPIN-BOSON MODEL
In this section, we turn to an explicit application of themethods developed above to the paradigmatic spin-boson
Figure 2. (a) Steady-state phonon emission spectrum for the Ohmicspin-boson model computed using the generalized transfer-tensormethod. The parameters used are α = 0 . , δt = 0 . (cid:15) − , k B T =0 . (cid:15) and ω c = 10 (cid:15) . The plot shows the spectra obtained bothwhen including three-time transfer tensors in the propagation (solidline), and the Markovian result obtained when these terms are ne-glected (dotted line). The time τ denotes the cutoff time used for thetransfer-tensor propagation, for τ = 0 . (cid:15) − the curve is offset by0.3 and for τ = 2 (cid:15) − by 0.6. (b) Relative entropy measure of non-Markovianity (c) Frobenius norm of the two-point transfer tensors.Notice the increase of the transfer tensor norm at long times. Thiseffect is well-known in the context of combined TEDOPA-TTMcomputations [26], where it signifies a back-reflection of excitationsemitted into the environment due to a finite-size effects. For thecase considered here the transfer-tensor method allows us to removefinite-size effects, by simply neglecting these transfer tensors. model [6]. First we apply the generalized transfer-tensormethod to study the steady-state emission spectra, and inves-tigate the effects of system-environment correlations presentin the steady state. We then make use of the fact that we arecomputing the process tensor, to operationally quantify thenon-Markovianity of the process [34], and study the dynami-cal build-up of system-environment correlations. A. Steady-state emission spectra
For concreteness we consider a spin-boson type model witha single two-level system interacting with a bosonic environ-ment. This system is characterized by the Hamiltonian H = ε σ x + σ z (cid:88) j g j (cid:16) a j + a † j (cid:17) + (cid:88) j ω j a † j a j (30)where σ x , σ y , σ z are Pauli spin operators, ε is the energysplitting of the two-level system, g j is a coupling coefficient, and a † j , a j are Bosonic creation and annihilation operators ofan environment mode with energy ω j . The effects of thebosonic environment on the spin dynamics is fully charac-terized by the spectral density J ( ω ) = αω c ( ω/ω c ) s e − ω/ω c ,where α is the coupling strength, ω c is the bath cutoff fre-quency and s is the Ohmicity. Two-, three- and multi-timedynamical maps can be efficiently simulated for this modelusing the TEMPO algorithm and its multi-time generaliza-tion [16, 36]. The reference state of the environment is takento be the thermal state with thermal energy k B T , where k B isBoltzmann’s constant.In many cases we are interested in studying second or-der correlation functions. This is the case when for exam-ple studying absorption and emission spectra of quantum sys-tems. In particular we can consider the steady-state emissionspectrum given by [37] S ( ω ) = Re (cid:20)(cid:90) ∞ dτ ( g (1) ( τ ) − g (1) ( ∞ )) e − iωτ (cid:21) , (31)defined in terms of the two-point correlation function g (1) ( τ ) = lim t →∞ (cid:104) σ + ( t + τ ) σ − ( t ) (cid:105) . The two-point cor-relation function is defined in terms of the raising and low-ering operators on the spin system σ ± = ( σ x ± iσ y ) . Wecan compute this function within our three-time measurementscenario, by implementing the raising and lowering operatorsin the place of measurement operators.In Fig. 2a we show the steady-state phonon emission spec-trum for the spin-boson model computed using the general-ized transfer-tensor approach. The spectrum has been cal-culated both with (black lines) and without (coloured lines)including the three-time transfer tensors propagating correla-tions across the implemented lowering operator. This corre-sponds to using an approximate set of transfer tensors, com-puted via T reg n,m,k = E n,m,k − T n,m ⊗ E m,k − n − (cid:88) j = m +1 T n,j E j,m,k . (32)This is equivalent to the quantum regression theorem ap-proach [6], in which it is assumed that no correlations arecarried across an implemented operation. In Fig. 2a we showthe computed spectra for three different memory cutoff times τ , at τ = 2 (cid:15) − the spectrum was converged, such that theblack line in this case gives the exact spectrum.In contrast to previous methods, the multi-time transfer-tensor approach applied here is able to fully account forsteady-state system-environment correlations. This is donewithout the need to consider inhomogeneous contributionsto the dynamics, as the three-time transfer-tensors makes itpossible to homogeneously propagate the state across an im-plemented operation. From the figure we observe that thesesteady-state correlations play an important role in the emis-sion spectrum, in particular we see that including the three-time transfer tensors tend to suppress the negative frequencypart of the spectrum. B. Quantifying non-Markovianity
Following the operational approach introduced by Pollocket al. [34], we define the sub-manifold of process tensor Choistates which can be expressed as product states, as Markovianprocesses. This class of processes represent the case where nocorrelations can be carried across an implemented operation.Having defined the Markov process, we can define a measureof non-Markovianity N . We do this using the quantum rela-tive entropy, which for the three-time process tensor gives N n,m, = tr (cid:2) Υ n,m, (cid:0) ln Υ n,m, − ln Υ markov n,m, (cid:1)(cid:3) . (33)Here the Markovian process is given by Υ markov n,m, = E n,m ⊗E m, ⊗ ρ . Measuring non-Markovianity as a relative entropyhas a clear operational significance as a hypothesis testingprocedure. That is how likely are we to confuse the Marko-vian process tensor for the non-Markovian process tensorgiven λ repetitions of the process. The probability of con-fusing the Markovian and Non-Markovian process tensors isgiven by P n,m, = exp ( − λ N n,m, ) , such that a larger rela-tive entropy gives a smaller likelihood of confusing the twoprocesses [34].In Fig. 2b we show the relative entropy computed for theohmic spin-boson model. Starting out from an initial productstate, we see how system-environment correlations are signif-icant at short times following an intermediate measurement.At longer times the system-environment correlations presentat the time of the intermediate measurement becomes negli-gible. This means that if we are interested in processes inwhich two implemented measurements do not occur to closetogether in time, then a Markovian description of the processis adequate. In the specific case of Fig. 2b, we see that anynon-Markovian effects are negligible when the second mea-surement is implemented at times greater than roughly (cid:15) − after the firts measurement. V. CONCLUSION
To conclude, we have presented a generalization of thetransfer-tensor formalism to include multi-time measure-ment scenarios. The generalization was naturally formulatedwithin the process tensor framework. We applied our methodto the problem of computing steady-state phonon emissionspectra for the spin-boson model. This problem has been ad-dressed before, the novelty of our computation is that we cantake full account of steady-state correlations without the needto compute inhomogeneous terms. Furthermore we couldstudy the build up of system-environment correlations duringthe dynamical evolution, by looking at the relative entropymeasure of non-Markovianity.In the future it would be interesting to investigate whetheran operational non-Markovianity measure could be related di-rectly to the transfer-tensor truncation error. Perhaps, thiswould make it possible to give notions such a memorystrength and process recoverability [38] a direct practical rele-vance as quantifiers of errors in the simulation of open quan-tum system dynamics. In addition, the methods presentedhere are likely to find many interesting applications. As an ex-ample we mention the study of emission and absorption spec-tra of organic and inorganic semiconductors, here one coulddirectly look at the effects of non-Markovianity, coherenceand temperature dependence.
ACKNOWLEDGMENTS
The authors would like to thank Jonatan B. Brask for help-ful discussions. MRJ was supported by the Independent Re-search Fund Denmark. [1] Weizhe Edward Liu, Jesper Levinsen, and Meera M. Parish,“Variational approach for impurity dynamics at finite tempera-ture,” Phys. Rev. Lett. , 205301 (2019).[2] Yuto Ashida, Tao Shi, Richard Schmidt, H. R. Sadeghpour,J. Ignacio Cirac, and Eugene Demler, “Quantum Rydberg cen-tral spin model,” Phys. Rev. Lett. , 183001 (2019).[3] Niclas Schlünzen, Jan-Philip Joost, and Michael Bonitz,“Achieving the scaling limit for nonequilibrium Green func-tions simulations,” Phys. Rev. Lett. , 076601 (2020).[4] Chloe Clear, Ross C. Schofield, Kyle D. Major, Jake Iles-Smith, Alex S. Clark, and Dara P. S. McCutcheon, “Phonon-induced optical dephasing in single organic molecules,” Phys.Rev. Lett. , 153602 (2020).[5] Alistair J. Brash, Jake Iles-Smith, Catherine L. Phillips, DaraP. S. McCutcheon, John O’Hara, Edmund Clarke, BenjaminRoyall, Luke R. Wilson, Jesper Mørk, Maurice S. Skolnick,A. Mark Fox, and Ahsan Nazir, “Light scattering from solid-state quantum emitters: Beyond the atomic picture,” Phys. Rev.Lett. , 167403 (2019).[6] H.-P. Breuer and F. Petruccione,
The Theory of Open QuantumSystems (Oxford University Press, 2002). [7] Inés de Vega and Daniel Alonso, “Dynamics of non-Markovianopen quantum systems,” Rev. Mod. Phys. , 015001 (2017).[8] Shaul Mukamel, “Femtosecond optical spectroscopy: A directlook at elementary chemical events,” Annu. Rev. Phys. Chem. , 647–681 (1990).[9] Joel Yuen-Zhou, Jacob J Krich, Ivan Kassal, Allan S John-son, and Alán Aspuru-Guzik, Ultrafast Spectroscopy , 2053-2563 (IOP Publishing, 2014).[10] Tobias Brixner, Jens Stenger, Harsha M. Vaswani, MinhaengCho, Robert E. Blankenship, and Graham R. Fleming, “Two-dimensional spectroscopy of electronic couplings in photosyn-thesis,” Nature , 625–628 (2005).[11] Crispin Gardiner and Peter Zoller,
Quantum Noise: A Hand-book of Markovian and Non-Markovian Quantum Stochas-tic Methods with Applications to Quantum Optics (Springer-Verlag Berlin Heidelberg, 2004).[12] Giacomo Guarnieri, Andrea Smirne, and Bassano Vacchini,“Quantum regression theorem and non-Markovianity of quan-tum dynamics,” Phys. Rev. A , 022110 (2014).[13] Dara P. S. McCutcheon, “Optical signatures of non-Markovianbehaviour in open quantum systems,” Phys. Rev. A , 022119 (2016).[14] Li Li, Michael J.W. Hall, and Howard M. Wiseman, “Conceptsof quantum non-Markovianity: A hierarchy,” Phys. Rep. ,1–51 (2018).[15] M. Cosacchi, M. Cygorek, F. Ungar, A. M. Barth, A. Vagov,and V. M. Axt, “Path-integral approach for nonequilibriummultitime correlation functions of open quantum systems cou-pled to markovian and non-markovian environments,” Phys.Rev. B , 125302 (2018).[16] Mathias R. Jørgensen and Felix A. Pollock, “Exploiting thecausal tensor network structure of quantum processes to ef-ficiently simulate non-Markovian path integrals,” Phys. Rev.Lett. , 240602 (2019).[17] Hsing-Ta Chen, Guy Cohen, and David R. Reichman, “Inch-worm Monte Carlo for exact non-adiabatic dynamics. i. theoryand algorithms,” J. Chem. Phys. , 054105 (2017).[18] Hsing-Ta Chen, Guy Cohen, and David R. Reichman, “Inch-worm Monte Carlo for exact non-adiabatic dynamics. ii. bench-marks and comparison with established methods,” J. Chem.Phys. , 054106 (2017).[19] Daniel Alonso and Inés de Vega, “Hierarchy of equations ofmultiple-time correlation functions,” Phys. Rev. A , 052108(2007).[20] Guy Cohen and Eran Rabani, “Memory effects in nonequi-librium quantum impurity models,” Phys. Rev. B , 075150(2011).[21] Bassano Vacchini, “Generalized master equations leading tocompletely positive dynamics,” Phys. Rev. Lett. , 230401(2016).[22] Qiang Shi and Eitan Geva, “A new approach to calculating thememory kernel of the generalized quantum master equationfor an arbitrary system-bath coupling,” J. Chem. Phys. ,12063–12076 (2003).[23] Javier Cerrillo and Jianshu Cao, “Non-Markovian dynamicalmaps: Numerical processing of open quantum trajectories,”Phys. Rev. Lett. , 110401 (2014).[24] Felix A. Pollock and Kavan Modi, “Tomographically recon-structed master equations for any open quantum dynamics,”Quantum , 76 (2018).[25] Andrei A. Golosov, Richard A. Friesner, and Philip Pechukas,“Efficient memory equation algorithm for reduced dynamics inspin-boson models,” J. Chem. Phys. , 138–146 (1999).[26] Robert Rosenbach, Javier Cerrillo, Susana F Huelga, Jian-shu Cao, and Martin B Plenio, “Efficient simulation of non-Markovian system-environment interaction,” New J. Phys. ,023035 (2016).[27] Maximillian Buser, Javier Cerillo, Gernot Schaller, and Jian-shu Cao, “Initial system-environment correlations via the trans-fer tensor method,” Phys. Rev. A , 062122 (2017).[28] Anton Ivanov and Heinz-Peter Breuer, “Extension of theNakajima-Zwanzig approach to multitime correlation func-tions of open systems,” Phys. Rev. A , 032113 (2015).[29] Felix A. Pollock, César Rodríguez-Rosario, Thomas Frauen-heim, Mauro Paternostro, and Kavan Modi, “Non-Markovianquantum processes: Complete framework and efficient charac-terization,” Phys. Rev. A , 012127 (2018).[30] Simon Milz, Fattah Sakuldee, Felix A. Pollock, and KavanModi, “Kolmogorov extension theorem for (quantum) causalmodelling and general probabilistic theories,” Quantum , 255(2020).[31] Fattah Sakuldee, Simon Milz, Felix A. Pollock, and KavanModi, “Non-markovian quantum control as coherent stochastictrajectories,” Journal of Physics A: Mathematical and Theoret-ical , 414014 (2018).[32] Ángel Rivas, Susana F. Huelga, and Martin B. Plenio, “Entan- glement and non-Markovianity of quantum evolutions,” Phys.Rev. Lett. , 050403 (2010).[33] Heinz-Peter Breuer, Elsi-Mari Laine, and Jyrki Piilo, “Mea-sure for the degree of non-Markovian behavior of quantum pro-cesses in open systems,” Phys. Rev. Lett. , 210401 (2009).[34] Felix A. Pollock, César Rodríguez-Rosario, Thomas Frauen-heim, Mauro Paternostro, and Kavan Modi, “OperationalMarkov condition for quantum processes,” Phys. Rev. Lett. , 040405 (2018).[35] Simon Milz, M. S. Kim, Felix A. Pollock, and Kavan Modi,“Completely positive divisibility does not mean Markovianity,”Phys. Rev. Lett. , 040401 (2019).[36] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W.Lovett, “Efficient non-Markovian quantum dynamics usingtime-evolving matrix product operators,” Nat. Commun. (2018).[37] Dara P. S. McCutcheon, “Optical signatures of non-Markovianbehavior in open quantum systems,” Phys. Rev. A93