Numerical evaluation of two-time correlation functions in open quantum systems with matrix product state methods: a comparison

We compare the efficiency of different matrix product state (MPS) based methods for the calculation of two-time correlation functions in open quantum systems. The methods are the purification approach [1] and two approaches [2,3] based on the Monte-Carlo wave function (MCWF) sampling of stochastic quantum trajectories using MPS techniques. We consider a XXZ spin chain either exposed to dephasing noise or to a dissipative local spin flip. We find that the preference for one of the approaches in terms of numerical efficiency depends strongly on the specific form of dissipation.


Introduction
The investigation of open quantum many-body systems has been a very active field of research over the past decades. One motivation is the understanding of destructive effects of environments on quantum processes which are used in quantum technologies, e.g in quantum computing or quantum communication. More recently, another point of view has been taken.
Environments are particularly tailored in order to stablize and control quantum many-body states [4][5][6][7][8]. However, it has been shown that the dynamic properties of such states can have very distinct behaviour from their Hamiltonian counterparts [9]. In particular the response of open quantum systems to external perturbations very different to the Hamiltonian evolution.
Quantities which are of particular importance in this respect are two-time correlation functions B(t 2 )A(t 1 ) . Here A and B are operators, t 1 and t 2 are two different times, and . . . = tr(ρ . . . ) is the expectation value over the density matrix ρ of a given system. In isolated systems, such two-time correlation functions are powerful tools to give information on the response of the system to a small perturbation. Many experimental techniques are based on such processes and the observation of the subsequent response is described by these two-time functions. Examples include neutron scattering [10], ARPES [11], conductivity and magnetization measurements in solids [12] and spectroscopic measurement as radio-frequency [13], Raman, Bragg [14] or modulation spectroscopy [15] in the field of quantum gases. Two-time correlations have been studied extensively in isolated many-body quantum systems both in and out of equilibrium. However, in many-body quantum systems coupled to environments their determination is very challenging and only few studies are available mostly using approximate approaches or small systems [9,[16][17][18][19][20]. In particular, a change of the behaviour of dynamic correlations has been demonstrated in the presence of dissipation (see for example [21][22][23][24][25][26]).
In this work we present a comprehensive study on the application of matrix product state (MPS) algorithms to the determination of two-time correlation functions in open systems. We compare an extension of the purification approach [27,28] to two-time correlations and two different stochastic approaches based on the unraveling of the quantum evolution [3,[29][30][31]. The first approach has been proposed by Breuer et al. [2] and the second approach by Mølmer et al. [3]. The comparison is performed using an XXZ spin model with two different couplings to the environment. The first coupling is a dephasing noise applied globally to the system and the second a local loss of magnetization. In section 2 we describe the models used. In section 3 we give an introduction to matrix product state based methods for open quantum systems. In section 4, we introduce the three different methods in order to calculate two-time correlation functions and in section 5 we present a comprehensive comparison between the methods.

Model
We consider spin-1 /2 chains, with coupling between spins on adjacent sites, as a paradigmatic model for interacting one-dimensional systems. In addition, the system is subjected to an environment which introduces dissipative processes to the system dynamics. Assuming that retroactive influences of earlier dissipation effects on the current dynamics can be neglected, the system dynamics can be described by the Lindblad master equation [31] ∂ ∂t with the superoperator L and the system density matrix ρ. The first term on the right hand side represents the unitary contribution generated by the Hamiltonian. Here we consider the XXZ spin-1 /2 Hamiltonian describing a chain of L spins, where J x and J z are exchange couplings according to different spin directions and S α l is the spin operator in direction α at site l. In equilibrium this model is well understood and exhibits in the ground state, three phases for different ratios of the interaction strengths [32,33]: For −1 ≤ J z /J x ≤ 1 a gapless Tomonaga-Luttinger liquid is formed, whereas J z /J x < −1 and J z /J x > 1 present gapped phases showing ferromagnetic and antiferromagnetic nature, respectively. The second term on the right hand side of Eq. 1 represents dissipative noise. We use the Lindblad form of the dissipator which is given by Here Γ l is the effective dissipation strength corresponding to the Lindblad jump operator L l . For the comparison of methods presented here, two different types of dissipation are considered. As will be shown, each of them have different impacts on various aspects of the performance of the methods. First we study systems exposed to bulk dephasing (see Fig. 1(a)), where the jump operators L l are given by the set of local S z l operators controlled by the dissipation strength Γ l = γ Local demagnetization < l a t e x i t s h a 1 _ b a s e 6 4 = " o q H T t w / W e E H k J P X h 0 y U u 6 V M n M c M = " > A A A C 9 3 i c j V L L T t t A F D 2 Y U h K g J c C y G 4 s I q a v I T t M A O 1 S g Y l M p S E 2 I l K B q 7 E y M G 7 8 0 d p D S i H / o t m z Y o W 7 7 O f x E v 6 C L n p k 6 l b q g 7 V j j e + f c c 8 + 9 8 / C y K M w L x 3 l Y s p a f r D x d r V T X 1 j e e P d + s b W 3 3 8 n S q f N n 1 0 4 s 7 t x 5 1 e s + H S P 2 / V j 9 6 U D 6 O C F 9 j F S 9 7 + P o 5 w h g 7 7 8 P E R n / E F t 9 b M u r P u r a + / q N Z S m b O D P 4 b 1 7 S d g c J p N < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o q H T t w / W e E H k J P X h 0 y U u 6 V M n M c M = " > A A A C 9 3 i c j V L L T t t A F D 2 Y U h K g J c C y G 4 s I q a v I T t M A O 1 S g Y l M p S E 2 I l K B q 7 E y M G 7 8 0 d p D S i H / o t m z Y o W 7 7 O f x E v 6 C L n p k 6 l b q g 7 V j j e + f c c 8 + 9 8 / C y K M w L x 3 l Y s p a f r D x d r V T X 1 j e e P d + s b W 3 3 8 n S q f N n 1 0 4 s 7 t x 5 1 e s + H S P 2 / V j 9 6 U D 6 O C F 9 j F S 9 7 + P o 5 w h g 7 7 8 P E R n / E F t 9 b M u r P u r a + / q N Z S m b O D P 4 b 1 7 S d g c J p N < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o q H T t w / W e E H k J P X h 0 y U u 6 V M n M c M = " > A A A C 9 3 i c j V L L T t t A F D 2 Y U h K g J c C y G 4 s I q a v I T t M A O 1 S g Y l M p S E 2 I l K B q 7 E y M G 7 8 0 d p D S i H / o t m z Y o W 7 7 O f x E v 6 C L n p k 6 l b q g 7 V j j e + f c c 8 + 9 8 / C y K M w L x 3 l Y s p a f r D x d r V T X 1 j e e P d + s b W 3 3 8 n S q f N n 1 0 4 s 7 t x 5 1 e s + H S P 2 / V j 9 6 U D 6 O C F 9 j F S 9 7 + P o 5 w h g 7 7 8 P E R n / E F t 9 b M u r P u r a + / q N Z S m b O D P 4 b 1 7 S d g c J p N < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " o q H T t w / W e E H k J P X h 0 y U u 6 V M n M c M = " > A A A C 9 3 i c j V L L T t t A F D 2 Y U h K g J c C y G 4 s I q a v I T t M A O 1 S g Y l M p S E 2 I l K B q 7 E y M G 7 8 0 d p D S i H / o t m z Y o W 7 7 O f x E v 6 C L n p k 6 l b q g 7 V j j e + f c c 8 + 9 8 / C y K M w L x 3 l Y s p a f r D x d r V T X 1 j e e P d + s b W 3 3 8 n S q f N n 1 0 Here the jump operators S z l are Hermitian, which implies that the infinite temperature state is a steady state of the model. One can show that this is the unique steady state. The dissipative dynamics arising in this and related models has been studied previously and interesting critical dynamics and aging dynamics have been pointed out to occur [1,18,[34][35][36].
Furthermore, we investigate the case in which the coupling to the environment results in a local defect that is represented by the single Lindblad operator S − c which only acts on the central site as shown in Fig. 1(b). Here c is the index of the central site for a chain with an odd number of sites and the center-left site otherwise. The dissipation strength is Γ c = γ and Γ l =c = 0 The steady state of this system is the ferromagnetic state with all spins pointing down. Using the mapping to interacting fermions or hard core bosons, the jump operator corresponds to a local particle loss process. Similar 'lossy' defects have been studied in a variety of models previously and interesting transport effects and metastable states have been identified [37][38][39][40][41][42][43][44][45][46].

Matrix product state formalism for closed systems
The description of quantum states in MPS form has become a standard method for the simulation of one-dimensional many-body quantum systems. It has been used for a wide range of models, since it is very efficient and well-controlled approximation [47]. In this section we describe the basics and key concepts of this technique.

MPS representation of quantum many-body states
The idea relies on the approximate representation of the quantum many-body wave function for a one-dimensional lattice system of L sites as a set of local tensors/matrices. In order to achieve this, a singular value decomposition (SVD) is applied to the amplitude matrix of a bipartite system whose subparts A and B are connected by a bond between site l and l + 1 Here |m A/B are basis states in subsystem A/B and U , s, V are obtained by the singular value decomposition of the amplitude matrix ψ m,n . Iterating this procedure on all bonds then yields the expression of the quantum state by single site tensors where we use | σ = |σ 1 σ 2 . . . σ L and σ l labels the local basis states at the site l ∈ {1, 2, . . . , L}.
The von Neumann entropy is a measure for the entanglement between the two subsystems connected by the considered bond. If the entanglement is not too strong, this corresponds to a sufficiently fast decay of the descendingly sorted squared singular values s 2 a l . In this case, the dimension of the matrices for the representation of the state can be cut at a maximal value D, resulting in a compressed state which is a very good approximation of the exact state. A weak von Neumann entanglement is for example found for ground states of shortranged one-dimensional Hamiltonians [47]. As the value D limits the extend of the indices {a 1 , a 2 , . . . , a L−1 }, connecting two adjacent site tensors, this parameter is known as the bond dimension. The approximation is controlled by the sum of the discarded squared singular values, the so-called truncation weight and reduces the computational complexity from exponential to polynomial. In the following we will use the established graphical notation for tensor networks [47] where a tensor is portrayed by a shape (here circle) with sticking out lines representing the indices. Connecting two lines depicts a tensor contraction with regard to this pair of indices. Using this notation, a MPS is depicted by

Time-dependent matrix product state algorithm
Time-dependent matrix product state (tMPS) algorithms are well-established tools for the efficient computation of the dynamics of closed many-body quantum systems [47][48][49][50]. In particular, it is possible to calculate the time evolution of a MPS for the considered spin-1 /2 model with γ = 0, using a Suzuki-Trotter decomposition [51,52] of the time evolution operator for small time steps ∆t. Here we use the second order Suzuki-Trotter decomposition given by Here H odd (H even ) only contains parts of the Hamiltonian covering the odd (even) numbered bonds of the lattice, so that all contributing terms are commuting, but [H odd , H even ] = 0. In the diagram formalism, this corresponds to a successive application of two-site gate operations and compressions where the application order is indicated by the dotted arrow and the different colors distinguish gates and subsequent compressions of H odd and H even . By contracting these gates with the MPS tensors, the bond dimension also increases from D to d 2 D, so that a subsequent compression via SVD is necessary. A large reduction of the computational effort can be achieved by including conservation laws. In the present case, the total magnetization M tot = j S z j is preserved by the closed system evolution under the Hamiltonian H XXZ . This leads to a block-diagonal form of matrices, and thus, tensor operations, such as SVDs or tensor contractions, can be performed much more efficiently.

Purification of the density matrix
One way to transfer the techniques from the previous section to dissipative systems with finite dissipation strength γ, is to rewrite the density matrix acting on the physical Hilbert space Figure 2: Dissipative evolution of the purified density matrix in MPS form by a single time step ∆t using a second order Suzuki-Trotter decomposition for the evolution operator represented by a consecutive application of four-site gates. The two colors mark the affiliation to one of the two parts of the Lindbladian, i.e. either L odd or L even , and the arrow indicates the order of application.
H phys as a state in a doubled space H phys ⊗ H phys [28,47,53] where |ρ ∈ H phys ⊗ H phys . Thus, the density matrix acting on the physical space becomes a pure state in the 'doubled' space, giving the procedure the name purification. The order of the resorting is in principle arbitrary. We use here the given resorting of the state, since this has the advantage that the time-evolution only acts on four 'neighbouring sites'. In this representation of the density matrix in the super space H phys ⊗ H phys , the Lindblad master equation (Eq. 1) is written as, with the new representation of the dissipator D in the super space, Here I is the identity matrix in the Hilbert space H phys and T denotes the transpose of an operator.

Representation of an initial state in the purified form
Often we are confronted with the situation that we start a Lindblad evolution from a pure state represented in the MPS formalism, for example this can be the ground state of a Hamiltonian obtained by the MPS ground state search. For this purpose we present the steps of purifying a state |ψ that is given in MPS form [Eq. 11]. Since the density matrix of this pure state is given by |ψ ψ|, first a copy of the state is created which represents the bra contribution. Then the tensors stemming from the ket and the bra corresponding to the same site are contracted, which results in a set of tensors with six indices In order to obtain single-site tensors for the ket and the bra part, the bond indices are combined, i.e. (a l−1 , a l−1 ) → λ l−1 and (a l , a l ) → λ l , before the tensor are separated into the respective single-site parts using a singular value decomposition with regard to the indices (λ l−1 , σ l ) × (σ l , λ l ). The new bond index, created by the SVD, is denoted by λ l , so that the MPS representation of a purified state reads This translates to the following diagrammatic expression It is important to notice that this approach affects the bond dimension drastically. Assuming the bond dimension of the state to purify is given by D, the combining of bond indices results in a dimension D 2 for the index set {λ l } and the SVD without compression causes an increase to dD 2 for the indices {λ l }. In order to keep the original specified bond dimension D, the state to purify needs to be compressed to a bond dimension √ D and the spectrum of SVD in the second step needs to be truncated after the D largest values. This poses a very strong constraint, so that particularly strongly entangled states are difficult to purify. Either the need for computational resources increases quadratically or the accuracy, measured by the truncation weight, becomes significantly worse.

Time-evolution of the purified density matrix
In analogy to the unitary closed system evolution (Eq. 13), it is also possible to split the Lindbladian into two contributions where L odd (L even ) covers the ket and bra sites connected by the odd (even) bonds. Again, the terms within L odd/even commute, but [L odd , L even ] = 0. Consequently, the dissipative evolution can also be approximated by the second order Suzuki-Trotter decomposition, which in this case results in a sequence of applications of four-site gates as shown in Fig. 2. This fact raises the challenge, that a gate application causes a stronger increase of the bond dimension The dimension of the individual indices is marked in blue color, revealing that the bond dimension of a purified state grows from D to up to d 4 D by applying a time evolution gate. This requires a stronger truncation to compress the bond dimension to the original size. Another important feature is that the use of quantum number conserving codes for the computation of the dissipative evolution is only possible if the Lindbladian is subject to a strong symmetry [54,55], i.e. not only the Hamilton operator but also every single jump operator of the model respects the conservation law. In the cases considered here, this only applies to D 1 (Eq.4), as the jump operator S − c in D 2 (Eq.5) does not conserve the total magnetization.

Calculating expectation values within the purification approach
Provided with the time evolved purified density matrix, we are left with the task to calculate expectation values of observables to extract information about the system. The trace relation for the expectation value of an operator A translates to a scalar product Here the first and second spin in |σ σ l denote the spin at site l for the ket and the bra part in the purified notation respectively. The state |1 is the purification of the unnormalized infinite temperature state [56]. Unfortunately, the possibility to be able to encode this state as a product state as shown in Eq. 23 gets lost, when using good quantum numbers, where the vector of the purified density matrix cannot be separated into contributions of single sites anymore. As a result, this state can exhibit a potentially large bond dimension. Another obstacle becomes apparent when exploiting good quantum numbers in the algorithm is that |1 spans over all symmetry blocks, so that a restriction to a certain quantum number sector makes a manual selection of basis states fulfilling this condition very cumbersome. Here it is helpful to realize, that the purification procedure is subject to a gauge degree of freedom. The measured expectation value is invariant under local unitary transformations of the states [47]. This fact can be used to construct a favorable representation of the bra space of the constituents of |1 , defined by To this end, the transformation U is chosen such, that the quantum number of the full initial state is equally distributed over local pairs of ket and bra states. For example, if we consider the M tot = 0 symmetry sector in the Lindblad model using the dissipator D 1 , one possible transformation is the Pauli matrix in x-direction U = σ x on the bra sites, so that we can rewrite the state |1 as a product of local spin pairs whereσ = −σ. This transformation automatically guarantees the restriction of the trace generating state |1 to the symmetry sector selected by the initial state. With this, the problems of selecting basis states respecting the quantum number conservation as well as the large bond dimension of an MPS representation of |1 are both solved. Also the initial state and the Lindbladian gates need to be adapted to be consistent with this gauge choice. With U = L l=1 (I ⊗ U ), the transformation relations are Supposing that the observable can be brought to an efficient tensor representation, we can then proceed to calculate the expectation value. In the exemplary and important case of a local observable the measurement is carried out by the tensor contractions outlined in Fig.3.
As explained before only the case with dissipator D 1 conserves the total magnetization. Thus we use this transformation in order to consider the good quantum numbers in the simulation. Applying this transformation to the Lindbladian leaves the Hamiltonian contribution unchanged and only changes the sign of the first part of the dissipator D 1 . The observables are calculated in the transformed basis (see Fig. 3).

Monte-Carlo wave function method
Instead of evolving the full density matrix, an alternative way is to compute the evolution of wave function trajectories in the original Hilbert space [3,29,30]. This is at the expense of the need for a sampling over many realizations due to the presence of stochastic processes, originating from the action of the environment on the system. This approach, known as the unravelling of the master equation, can be realized by piece-wise deterministic jump processes, where the deterministic time evolution of a state is interrupted by the application of jump operators. The deterministic evolution is performed with regard to an effective Hamiltonian where (b) apply the selected operator L l to the state and renormalize it 5. Iterate from 2 until the final time is reached.
Performing the average over many trajectories generated by this scheme ultimately yields an estimate of the density matrix that is accurate up to the first order in the time step. The density matrix can be approximated by a Monte-Carlo (MC) average over a finite number R of trajectory samples |ψ r (t) giving the technique name of Monte-Carlo wave function (MCWF) method. Similarly, the expectation values of observables can also be evaluated as where . . . denotes the MC average and the stochastic accuracy of the independent sampling approach is given by the time-dependent standard deviation of the mean  The random sampling of an application time in steps 2 and 3 of the algorithm is equivalent to drawing a waiting time τ until the next jump occurs according to the distribution This quantity is important for the convergence of the method as it can be used to determine a sufficiently small time step when simulating the non-unitary evolution under the effective Hamiltonian. It is important to choose a time step small enough so that the probability of more than one jump event happening in the span of a time step is negligible. For the system with dissipator D 1 the imaginary part of the effective Hamiltonian (Eq. 27) is equal to − γ 2 l S z l S z l = − γL 8 . Thus the waiting time distribution (Eq. 33) simplifies to The waiting time distribution is the probability of having a jump in the interval [t, t + τ ) which is independent of the initial time t here. The probability to have two or more jumps in this interval is P n>2 (τ ) = 1 − e −γLτ /4 − γLτ 4 e −γLτ /4 . In Fig. 4 we show the probability of two or more jumps taking place in one time step for the model with dissipator D 1 . When a good value for ∆t has been found for one parameter configuration, the dissipation strength and the system size are crucial for estimating a suitable time step for another set-up.
For the determination of a single trajectory applying the introduced tMPS algorithm from Eq. 13 offers a promising solution for computing the deterministic part of the evolution [57]. Using this procedure, we can access the dissipative time evolution of a spin-1 /2 chain under the action of the Lindbladian. It is noteworthy that here the use of symmetries to reduce the computational complexity only requires the conservation of the associated quantum numbers by the effective Hamiltonian. Consequently, in contrast to the purification approach, quantum number conserving codes can be used for the simulation of both disspators D 1 and D 2 . In Fig. 5 we present for the example of the dissipator D 1 the evolution of the local equal-time correlation function Here, we find that during the evolution a spreading around the mean develops so that the maximum entanglement in individual samples can be significantly larger than the mean value suggests. The relationship between entropy and the necessary bond dimension is highlighted further for a single trajectory created from the same random seed for different maximal bond dimensions D in Fig. 6. It is evident that calculations with a fixed maximum value for the bond dimension can only provide a sufficiently accurate representation of the entanglement up to a certain time.

Methods for the computation of two-time correlation functions in open quantum systems
We will discuss the different concepts for determining two-time correlation functions in open quantum systems with the methods introduced in Sec. 3, which will later form the basis for the comparison. More precisely, the goal is to compute correlation functions of the form where the two operators are applied at different times during the course of a dissipative time evolution.

Two-time correlations within the purification approach
The purification approach can be straight-forwardly extended in order to determine two time correlation functions [1]. After purifying the initial state, the reshaped density matrix is evolved until time t 1 , where the operator A is applied, followed by an evolution to t 2 , where an ordinary measurement of B is carried out Only a single time evolution needs to be calculated. Nevertheless, the fact of the evolution taking place in a doubled Hilbert space, can potentially result in the need for a large bond dimension.

Stochastic sampling: two approaches to two-time correlations
In the following we introduce two different approaches to evaluate two-time correlation functions in the framework of Monte-Carlo wave functions. The first step is to evolve an initially prepared state, which is in the case of a mixed state drawn according to the weights in the density matrix, up to the application time t 1 of the first operator. This is performed using the unravelling scheme of piece-wise deterministic jump processes. For each trajectory one defines the state |φ(t 1 ) ≡ A |ψ(t 1 ) . In principle it seems that we are left with the task of calculating the dissipative time-evolution of the two states |φ(t 1 ) and |ψ(t 1 ) up to time t 2 and then the expectation value with the operator B, i.e.
However, while this expression is well-defined for states of a closed system, the transfer to a stochastic sampling approached is more involved. In the next two sections we outline two methods which are well defined for the stochastic sampling approach. Subsequently, we compare the efficiency of both concepts.

Joint evolution of two states following Breuer et al. [2]
The first idea, developed by Breuer et al. [2], uses a doubling of the Hilbert space at time t 1 by introducing for each trajectory the vector As the matrixρ(t 1 ) = |Θ(t 1 ) Θ(t 1 )| again fulfills all properties of a physical density matrix, the task is reduced to recover the time dependence in Eq. 38. By defining a new Hamiltonian operator and new jump operators acting on the doubled spacẽ and using this in a Lindblad-type equation with HamiltonianH and jump operatorsL l , we arrive at Lindblad equations for all matrix blocks [2]. As a result, we can apply the same unravelling approach as before on the doubled space to compute two-time correlation functions. Since the operators do not couple the two subspaces, a separate time evolution of |ψ(t) and |φ(t) is possible, where only the application time and the selection of jump operators is determined based on the joint evolution (see Fig. 7). The algorithm to create a single trajectory can be condensed to the following steps: 1. Initialize the wave function |ψ(t 0 ) in the original Hilbert space and evolve it until t 1 using the introduced piece-wise deterministic process.

Evolve both states independently under the effective non-Hermitian Hamiltonian while
sampling jumps simultaneously according to the joint loss of norm of |Θ(t) .
Then, use the components of this state in order to measure the two-time correlation function by calculating the full overlap This method can be further extended to multi-time correlation functions by additional operator applications between the first and the final application time [2].
split trajectories joint jumps L j Figure 7: Sketch of the creation of one sample for the computation of two-time correlators following [2]. After the evolution up to t 1 the trajectory is copied, followed by a time span characterized by an independent deterministic non-unitary evolution interrupted by joint jump operator applications.

Separate evolution of four trajectories following Mølmer et al. [3]
An alternative strategy, established by Mølmer et al. [3], uses the quantum regression theorem [31] to represent the time dependence of two-time correlation functions in terms of the evolution of quantities, which only involve equal-time measurements. The proof relies on the quantum regression theorem which states that if the evolution of equal-time observables for the set of operators {B i } is described by a closed set of differential equations the two-time correlations are generated by the same kernel G as The idea of this approach relies on constructing four new states at time t 1 for each trajectory such that a combination of the expectation values of the operator applied at t 2 for each of these state gives the corresponding two-time correlation. To this end, the four states are defined by applying the first operator A [3], where for each trajectory |ψ(t 1 ) the state at t 1 is obtained by the usual unravelling scheme. Further, µ ± R and µ ± I normalize the respective states. These new states evolve independently from time t 1 to t 2 with the same unravelling scheme (see Fig. 8). By properly combining the expectation values of the second operator B for each state the two-time correlation functions can be recovered [3] Consequently, it is possible to access the two-time correlation functions by evolving the four states from Eq. 43 separately. In contrast to the MCWF technique of the doubled Hilbert independent evolution Figure 8: Sketch of the creation of one sample for the computation of two-time correlators following Ref. [3]. The initial dynamics up to time t 1 is followed by the splitting of the wave function into four trajectories from which the two-time expectation value at t 1 can be reconstructed. Then each sub-trajectory is evolved separately.
space, also the jump sampling procedure is completely independent so that it is possible to compute the four trajectories after t 1 in parallel. In addition to the average over the four states, the average over many trajectories from the initial time to t 1 and from time t 1 to t 2 needs to be taken. A disadvantage of this approach using four different state is, however, that it does not offer a straight-forward extension to multi-time correlation functions. There are two special cases, defined by the properties of the operators A and B of the correlation function. If the operator at time t 1 is Hermitian (A † = A) and commutes with the operator at t 2 ([A, B] = 0), it is possible to consider two instead of four trajectories and to calculate the two-time correlation as To ensure that the comparison between the two approaches is as fair as possible, we concentrate on such a special case for the two-time correlator in Eq. 46. The second case is more subtle and appears in the context of conserved quantum numbers. If the Lindbladian evolution allows the partition of the superoperator into symmetry blocks, and the operator A at t 1 couples different blocks, the resulting states |χ ± R,I cannot longer be addressed to a single block. Instead it is possible to evolve the parts of the different sectors separately as they are decoupled for all times t > t 1 . A closer look reveals that in this case the introduced scheme is equivalent to the approach by Breuer et al. from the last section.

Comparison of the different methods to determine two-time correlations functions in open quantum systems
In this section we compare the performance of the different methods to calculate two-time correlations using the spin model described in section 2. We focus on the two-time correlation function relating two applications of local S z -operators at sites separated by a distance d given by Recall, c is the central site for a chain with an odd number of sites and the center-left site otherwise. This correlation function has been proven to be essential to uncover interesting dynamical regimes displaying physical phenomena such as aging or hierarchical dynamics [1]. We start by the comparison of the two stochastic approaches for the dissipator D 1 in subsection 5.1. We use as an initial state a classical Neel state. We find that typically the Breuer et al. approach combined with the MPS methods performs better than the Mølmer et al. approach. We continue in comparing the better performing stochastic approach, the approach by Breuer et al., to the purification approach in subsection 5.2. For the considered situation the purification approach greatly outperforms the stochastic approach. Reasons of the excellent performance of the purification approach in this situation are the 'easy' initial state and the conservation of the magnetization by the Lindblad dynamics and most importantly, the low matrix dimension needed. In subsection 5.3 we turn to the local dissipator D 2 . We consider the situation where the initial state is the ground state of the Hamiltonian and the dissipation is switched on at time t = 0. We compare again the stochastic approach of Breuer et al. to the purification approach. The entangled initial state and the abscence of the conservation of magnetization constitutes an additional difficulty for the purification approach. We find that for most parameters considered the stochastic approach is much more suited to treat this situation efficiently.

Comparison of stochastic approaches for the dephasing noise D 1
In this section we compare the efficiency of the two stochastic aproaches in order to calculate two-time correlation functions.
We calculate the correlation function from Eq. 46 for the spin model with the dissipator D 1 . The choice for the correlation functions allows us to only use two states, i.e. |χ ± R , in the  First, we investigate the cost of generating a single two-time trajectory sample with MPS methods. To this end, the entanglement entropy is evaluated for all times t > t 1 . As it is directly linked to the needed bond dimension, it serves as an architecture-free measure of the required resources in terms of run time and memory consumption. In Fig. 9 we show the von Neumann entropy for the two branches in the Breuer approach, where |φ(t 1 ) = A |ψ(t 1 ) , as well as for the |χ ± R branches of the Mølmer procedure. As the results for the two Mølmer branches are indistinguishable in the plot, we only present data for |χ + R . Besides the MC average we also show the maximum measured entropy over all sampled trajectories and the 1σ interval of the measured standard deviation. In order to avoid biases introduced by the truncation we use an exact MPS representation for L = 14 sites. We see that the entropy increases over time and the statistical deviation increases in all cases. However, for both parameter sets presented it is evident that the two branches of the Breuer approach generate significantly less entropy than the Mølmer branches. Furthermore, a strong dependence on the model parameters exists. The von Neumann entropy is much smaller for large anisotropy parameters J z /J x .
Another important factor is the convergence of the Monte-Carlo averages with respect to the number of samples. We present the scaling of the standard deviation of the mean with the number of samples at different time points in Fig. 10. As the sampling is statistically independent we observe an inverse square root scaling in all cases. However, this is more than one order of magnitude smaller for the joint evolution suggested by Breuer et al. than for 10 0 10 1

Comparison of Monte-Carlo wave function and purification method for D 1
After having identified the approach by Breuer et al. as superior compared to Mølmer et al. for the considered situation, we compare this approach to the purification method. Since the exact increase of the bond dimension is unknown in both cases, the trade-off between the larger Hilbert space and the cost of sampling many trajectories needs to be evaluated carefully.
We compare the accuracy of the two-time correlation function C d (t 2 , t 1 ) normalized to the value at t 1 during the dissipative evolution of a spin-1 /2 chain initially prepared in the Neel state with bulk dephasing, i.e. using the dissipator D 1 . We previously performed such calculations in Ref. [1] and analyzed there the time step required to obtain a good accuracy. In the following we use these time steps for our comparison. Fig. 11 shows the results obtained by the purification method with different values for the bond dimension. This method converges faster in terms of the bond dimension for larger values of the dissipation strength and smaller system sizes. Since it is not clear at which time the greatest inaccuracy exists, we use the maximum deviation of the normalized two-time correlations from the results of the largest achievable bond dimension D max over the entire time as a measure for the error The two-time correlation C 1 (t 2 , t 1 ) which is calculated using bond dimension D is denoted by C 1 (t 2 , t 1 , D). The results for the stochastic sampling approach presented in Fig. 12 point in a similar direction, i.e. stronger dissipation and smaller systems require a smaller bond dimension to yield reliable results. Surprisingly, the convergence of the stochastic approach with increasing the bond dimension is generally much slower compared to the purification method. In addition, the accuracy is also influenced by the number of stochastic samples taken. To capture this effect, we first define the maximum value of the standard deviation of the mean in time as and then the total error for the MCWF approach as the maximum of ∆ max (D) and σ max (D). As can be seen in Fig. 13, depending on the number of samples used, the sampling error   surpasses the error caused by the MPS truncation when using a sufficiently large bond dimension.
Another important fact to consider is the dependence of the choice for the time step on the dissipation strength and the system size in the MCWF approach as mentioned in the discussion of Fig. 4. Fig. 14 confirms that ∆tJ x / = 0.02 is a suitable time step for the dissipation strength γ/J x = 1. Using the extrapolation indicated in Fig. 4 we estimate ∆tJ x / ≈ 0.002 as a potential step size for γ/J x = 10. However, the convergence analysis reveals that this is still not sufficient to obtain high quality simulation results. We find that larger dissipation strengths require a very small time step for MCWF simulations in this cases, which causes substantially longer runtimes and additionally increases the error caused by the MPS truncations, as more gate application with subsequent singular value decompositions are needed to reach a certain simulation time.
With the introduced benchmark strategy, i.e. using Error(D) = ∆ max (D) : purification max {∆ max (D), σ max (D)} : MCWF (49) as an error measure for the two approaches, it is now possible to compare them quantitatively. The dependence of the error on the bond dimension is presented in Fig. 15 for two different system sizes and dissipation strengths. While a weak dependence on the system size is noticeable, the dissipation strength turns out to be the decisive parameter. For the purification calculations the convergence regarding the bond dimension is much faster for strong dissipation. The accuracy of the stochastic sampling is relatively quickly dominated by the statistic error, so that the behaviour with the bond dimension appears to be almost constant. Consequently, there are crossing points above which the accuracy of the purification method is better for the same bond dimension.
To put these results into the context of realistic numerical resources, we show in Fig. 16 the required run time for a certain error on a computer cluster consisting of machines with 2.6 GHz clock frequency and sufficient memory resources. In case of strong dissipation, the generation of a single trajectory already takes longer than the evolution using the purification of the reshaped density matrix due to the necessity of a much smaller time step. For weaker dissipation strengths the run time of a single trajectory becomes comparable to the purification approach. However, the total run time of the stochastic sampling computation, requiring a sample size of the order 10 4 trajectories, is orders of magnitude larger than the full evolution, even with a massive parallelization of the sampling process. This means that the purification approach in these cases is strongly favourable over the stochastic approach.

Comparison of purification and stochastic approach for local dissipation D 2
To demonstrate the strong influence of the model and the initial state on the method choice, we continue by supplementing the findings of the last chapter with investigating results of another set-up. For this purpose, we turn to the dissipator D 2 , which only contains one jump operator given by the spin lowering operator S − c acting on the central site of a system with an odd number of spins. As this jump operator violates the conservation of the total magnetization during the Lindbladian evolution, it is not possible to exploit symmetry properties in the evolution of the purified state. Nevertheless, the non-unitary evolution in the deterministic part of the stochastic sampling conserves the total magnetization and only the jump applications switch between different symmetry blocks, so that the MCWF evolution can be calculated using conserved quantum numbers.
While the initial state was a product state in the previous chapter, the system is here initially prepared in the ground state of the equilibrium model accessed with the density matrix renormalization group (DMRG) ground state search algorithm [47,58]. We have chosen the ground state of an anisotropy of J z /J x = 0.5, which is located in the gapless Tomonaga-Luttinger liquid phase of the Hamiltonian and requires a sizable bond dimension for its representation in MPS form. As a result, a strong truncation is necessary when reshaping the corresponding density matrix to a purified state in the doubled Hilbert space as described  Figure 17: Deviation of ground state energy after purification compared to the value obtained by DMRG in the original Hilbert space. The ground state has been calculated for J z /J x = 0.5.
in Sec. 3.2.1. In Fig. 17 we show the bond dimension dependence of the deviation of the ground state energy calculated after the reshaping process from the original value obtained by DMRG. For small system sizes with less than 20 sites, the ground state energy can be reproduced after the purification step with a medium sized bond dimension with less than 10 3 states taken into account. However, even for slightly larger systems with up to 50 spins bond dimensions of several thousand states are needed to achieve an accuracy of the ground state energy of only 10 −3 . In this case the purification step alone takes more than three days of run time.
To estimate the numerical effort and the influence of inaccuracy in the representation of the ground state on time evolution for the purification method we compare the equal-time correlation functions for different values of the bond dimension in Fig. 18. As the equal-time correlations represent the initial condition for the two-time correlation functions at time t 1 , the accuracy of their calculation is crucial in order to obtain the two-time correlation functions. The direct comparison in Fig. 18 shows that the purification method reaches a comparable accuracy to the MCWF approach for D ≥ 600. Even for these bond dimensions sizable deviations occur of the order of 10 −3 . Looking at the associated run times, as summarized in Tab. 1, shows that the simulation based on purification takes about ten times as long for such a large bond dimension as compared to a parallel implementation of the stochastic sampling.
Based on the substantially larger run time, the purification becomes very inefficient and the MCWF scheme is preferentiable already to reach the time t 1 . To conclude the analysis, it remains to be demonstrated that the two-time correlation functions are accessible by the MCWF scheme. For this purpose we show in Fig. 19 the two-time correlation functions for the system sizes L = 25 and L = 29 which would be very inaccurate and time consuming to reach by the purification method.     Fig. 18. The MCWF simulation was executed in parallel on 10 cores. Parameters are the same as in Fig. 18.

Conclusion
To conclude we have presented a comparison of three different MPS based methods for the calculation of two-time functions in open quantum systems. This comprises the purification approach and two different approaches based on the stochastic unravelling of the Lindblad dynamics. First we compared the two stochastic approaches in the situation of an XXZ spin chain subjected to a dephasing noise starting initially in the classical Neel state. In this situation we find a clear preference for the stochastic approach suggested by Breuer et al. over the approach suggested by Mømler et al.. This is due to the better convergence of the trajectories used in the approach by Breuer et al.. However, the purification approach is even much more efficient for the considered situation. Additionally, we considered the dynamics of a XXZ spin chain subjected to a local application of the jump operator S − c starting from a Tomonaga-Luttinger liquid. This changes drastically the efficiency of the methods and the stochastic approach becomes more efficient than the purification approach. There are several reasons for this. First, already the representation of the Tomonaga-Luttinger liquid in the purified form is resource demanding. Secondly, in the following time-evolution the conservation of the magnetization is not fulfilled anymore.
We therefore expect that the purification approach is valuable if the initial state is easily represented within the purified space. Further, a strong symmetry of the Linbladian enabling the use of conserved quantities is of advantage. In comparison the stochastic wave function approach is well suited also to represent difficult initial states and the following Lindblad evolution. However, the presence of many jumps, as for the case of strong dephasing noise, calls for a very low time-step which makes the trajectory approach less efficient.