{'fraction_non_alphanumeric': 0.042300396046194295, 'fraction_numerical': 0.010699903484540887, 'mean_word_length': 4.277972949235903, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Systematic Gorini, Kossakowski, Sudarshan and Lindblad Equation for Open Quantum Systems

In complex, open quantum systems, with many degrees of freedom, it is difficult to start a Markovian dynamics, let alone prepare a factorized system-environment state. Rather, the Markovian dynamics is an intertwining map of a completely positive (CP) non-Markovian dynamics. The Markovian master equation is usually not CP as such, and should be left alone if it works adequately well. Here we study how the intertwining Markovian dynamics can still be well approximated by a CP map. A coherent quantum dynamics can be systematically approximated by the CP, Geometric Arithmetic Master Equation (GAME), but this is not the case for an incoherent dynamics. This dichotomy is responsible for disagreeing claims if a CP Markovian master equation cannot have systematic accuracy. As an example, we show how the dynamical decoupling of a qubit breaks the perturbative, intertwining map, while fixing the accuracy of the GAME to near-exactness.

Understanding the dynamics of open quantum systems is critical for quantum technologies, for, the openness to the environment, or to the bath, leads to a dephasing of quantum states. Models of open quantum dynamics sometimes display complex problems, ranging from negative probabilities of observables [1][2][3], instabilities of the reduced states versus time [4], to quantum phase transitions driven by the environmental coupling [5][6][7]. These problems generally do not happen if the reduced system dynamics is governed by a Gorini, Kossakowski, Sudarshan and Lindblad (GKSL) master equation (ME) [8,9]. The equation is rigorous only in the singular limit, where the bath correlation function (BCF) is the delta; also in the secular approximation (SA), which applies to an arbitrary BCF [10].
The defining feature of the GKSL ME is that it preserves CP of the reduced dynamics, in compliance with the axiomatic quantum mechanics. Here, we define a ME to be Markovian, if the coefficients of the equation are independent of the initial time where the systemenvironment state is factorized. Otherwise, the ME will be non-Markovian. A Markovian ME will be CP, iff it has the GKSL form [8][9][10]. The most important benefit of CP is that, by introducing a great mathematical simplification [11], CP becomes the "factory floor" in open quantum systems theory.
The underlying debate ranges from universality [15,31] to fallacy [32] of the ME with enforced CP. The nonnecessity of the GKSL form of a Markovian ME was first * dragomir.davidovic@physics.gatech.edu articulated by Whitney [33], but the debate is still ongoing [14, 15, 21-25, 32, 34-39]. The issue is, in a nutshell, that the bath always has a finite correlation time. So quantum dynamics always generates some entanglement between the system and the bath. Then, if the correlation time is much shorter than the time of the state of the system to change, in the interaction picture, the ME will not be Markovian initially, but likely become so after some time. Thus, the Markov property can only be asymptotic. Now, a Markovian master equation will be CP iff the initial system-environment state is factorized (or has a zero quantum discord [40]). Since the initiation of the Markovian dynamics into the classical state is thus impossible, there is no requirement for it to be CP.
Let us define a few terms, before we continue. In the following, we split the reduced density matrix of the system between the diagonal and the off-diagonal component, e.g., the populations and the coherences, respectively. To measure them, we utilize the respective trace norms; to measure their absolute errors, we use the respective trace distances to the exact populations and the exact coherences. The relative errors will be the ratios of the trace distances and the respective trace norms. In subsequent text, we will use the symbols ρ and to represent the state in the Schrödinger and the interaction picture, respectively.
If the relative errors are at a specified order in the coupling rate to the bath, for any quantum state, we will refer to the ME that solves them as universal. The attention to the relative error is important, because, depending on the state, the absolute and the relative errors can have different scaling in the coupling rate (see Sec. IV).
We define a ME to be systematic, if the relative errors of both the populations and the coherences are proportional to the coupling rate, for any coherent state, but not necessarily incoherent states. In the limit when the coupling rate goes to zero, the incoherent states span the range of quantum states of measure zero. Thus, a systematic ME can determine the coherent dynamics with arbitrary accuracy, by reducing the coupling rate to the environment. On the other hand, a universal ME can have that accuracy by increasing the order of the perturbative expansion, without reducing the coupling rate.
Let us restart the discussion with the question if there can be a systematic GKSL master equation? We first announced and subsequently published such an equation in Ref. [14]. In that work we derive a GKSL master equation by approximating the RE utilizing the geometric-arithmetic mean approximation, and name it the Geometric-Arithmetic Master Equation. The essence of the GAME, to be articulated in Sec. II, is that it advances the domain of the SA to a much wider range of system frequencies, without making any assumption about the distribution of the frequencies.
In Ref. [14] we determined the accuracy of the GAME by numerical simulation of the relaxation process from a macroscopic superposition in a Heisenberg ferromagnetic spin chain. Namely, for the initial magnetization direction perpendicular to the easy axis of the chain, the accuracy of the ensuing dynamics was determined by the trace-distance between the solutions of the RE and the GAME. Our key result was that the accuracy was proportional to the coupling rate to the bath [14], which is the same as the theoretical accuracy of the RE. It is important that our claim was based on numerical simulations, making the claim undisputable, although at expense of somewhat narrower scope as we shall learn here.
Shortly after the announcement of the GAME, a very similar master equation named the Universal Lindblad Equation (ULE) was announced by Nathan and Rudner [15,31], who claim universal accuracy based on a rigorous derivation and a rigorous error bound. The GAME and the ULE have identical dissipators, but differ in the unitary component of quantum dynamics known as the Lamb shift. The authors of the ULE make several claims that require more scrutiny: 1, that the ULE is derived by perturbative expansion in the coupling rate to the bath, and has perturbative accuracy; 2, that it is valid solely on properties of the bath and system-bath coupling strength, while independent of any energy level spacings in the system; and 3, that there exists an infinite family of distinct Markovian master equations valid on an equivalent level of the approximation.
Regarding the claim 1, a preprint announced recently by Tupkary et al. [32] analyzes the properties of the steady state solutions of the RE. They show that if any GKSL form is enforced onto the RE, the relative errors in the coherences become at zeroth order in coupling to the bath, contradicting the claims of universality [15,31]. They propose that the solution to all these issues must involve the forth order cumulant. We will show here that is not the case (see Sec. III, Eqs. 62 and 64).
Let us restate, from the outset, as in [32], that any claim of a universal GKSL equation, defined by the property that it solves both the populations and the coherences of reduced states with perturbative accuracy in the coupling rate, is a physical impossibility. For, we will show that the exact dynamics, of the time convolution-less master equation (TCL) [41], at infinite order in the coupling, does not have the GKSL form. Since any attempt to enforce the GKSL form makes the coherences non-systematic [32], the claim of universality is not in accordance with the facts. Claim 2 concerns the region of validity of the ULE as independent of the system frequencies. We will explicitly demonstrate the dependence of the error of the GAME on the oscillation frequencies of the system, and illuminate the rationale in Sec. II. The unfortunate flaw in the claim 2 is that it is based on a rigorous but weak error bound, terminating rather than emphasizing the frequency dependence (see Eq. 38 and related discussion for a proper error estimate).
Claim 3 proposes an infinite family of Markovian master equations with equivalent level of accuracy, including, but not limited to, the RE, the ULE, and the coarsegrained master equation (CGME) [25]. Regarding the latter, its accuracy is proportional to the square root in the coupling rate, in contrast to the RE that has linear accuracy. Thus, the RE and the CGME are not on the equivalent level. Similarly, as we show in this paper, when applied to incoherent energy transport, neither the ULE nor the GAME have accuracies at that level.
As an example of the usefulness of the GKSL enforced ME, in this paper we study the dynamical decoupling of a quantum system immersed in an Ohmic environment at zero temperature. We find that the (GKSL-enforced) GAME approaches exactness when the dephasing is suppressed, while the perturbative RE becomes unreliable.
The paper is organized as follows. In Sec. I we introduce the mathematical framework and the notation. In Sec. I A we determine the unique Lamb shift and the unique relaxation tensor in an arbitrary trace preserving and Hermiticity preserving ME. We rederive the GAME in Sec. II in wider regime of open quantum dynamics than in our prior work, discuss the "grey zone" of the equation, and emphasize the context of an advanced secular approximation. This is followed by Sec. III, where we identify non-CP Markovian dynamics as an intertwining map of a CP Non-Markovian dynamics. An unexpected but great benefit to the nonpositivity of the TCL equation at finite order in the coupling rate will be illuminated. The main results are presented in Sec. IV, where we study the systematic accuracy and the absence of the universal accuracy of the GAME or the ULE. Lastly, in Sec. V we discuss the merits of the GAME in dynamical decoupling.

I. PERTURBATIVE MASTER EQUATIONS
We start by summarizing the derivation of the RE adopting the notation in Ref. [41]. Consider a quantum system S weakly coupled to an environment B described by the Hamiltonian where H S and H B are the Hamiltonians of the system and the environment, respectively, and H I is the interaction Hamiltonian between the two. In the Schrödinger picture, H S can be time-dependent, but H I and H B are assumed to be time-independent. The von Neumann equation of motion for the density matrix of the total system is given in the interaction picture as This equation can be formally solved as where T is the chronological ordering operator. If we assume factorized initial condition of the system plus the environment, take the partial trace of Eq. 3 over the environment, and assume that the odd powers of H I are traceless in the environment, we arrive at the exact expansion in cumulants, which corresponds to [41,Eq. (9.65)] after tracing the environment. For the initial reduced state of the environment B (0), which from now on we write as B , we assume that is invertible as a linear map, the time convolutionless ME is derived by taking the derivative of Eq. 4, followed by expressing S (0) in terms of S (t) by the inverse map. From now one we adopt the notation TCL and TCL t , for the exact time convolutionless equation in the long time limit and at finite time, respectively. TCL t equation can be derived perturbatively in a systematic way, by utilizing the technique of time ordered cumulants, originally derived by van Kampen [42,43], or by the projection operator technique from the Nakajima-Zwanzig equation [44,45]. We use the notation TCL2n and TCL2n t if the equation is calculated at 2n-th order in the interaction with the bath. The RE is the time convolutionless equation TCL2. From now one we adopt the notation RE and RE t , for the long time limit and the finite time, respectively. The ME discussed so far are universal.
Let us begin with the form of the RE in Ref. [41,Eq. (3.117 which can be obtained from Eq. 4 by taking the time derivative, dropping the quartic commutator and substituting S (0) on the RHS with S (t).
We assume for the interaction Hamiltonian where A α and B α are Hermitian operators of the system and environment, respectively. In the interaction picture, we express A α (t) in terms of the Fourier transform, If H S is time-independent, with discrete eigenenergies E n and projector onto the eigenspace with eigenenergy E n denoted by Π n , A α (ω) will be determined from the eigenoperator decomposition [41], as If the system Hamiltonian is time-dependent, the frequencies will be quasi-continuous and the sum in Eq. 7 is to be replaced by the appropriate integral, where In either case, after some algebra, the RE can be written in the form [41,Eq. (3.132)], e.g., Note that this is not precisely the same as [41,Eq. (3.132)], because we did not extend the integration limit in Eq. 5 from t to ∞. Hence the time dependence in the spectral density matrix Here C αβ (s) = tr B [ B B α (s)B β (0)] is the bath correlation function (BCF). We also define the timed BCF as C αβ (s, t) = C αβ (s)Θ(s − t), so that the spectral density matrix is the half-sided Fourier transform of the timed BCF.
The RHS of Eq. 11 can be written down in a more useful form as where we define In comparison to the RE, the TCL master equation is very cumbersome, so we will leave all the details to appendix A.

A. Lamb Shift and Relaxation Tensor in Weakly Coupled Open Quantum Systems
In this section we extract the unique traceless Hamiltonian and the unique relaxation tensor (also known as the Kossakowski matrix) from an arbitrary trace preserving and Hermiticity preserving ME. In particular, for the RE in the interaction picture, the unique traceless Hamiltonian, or the Lamb shift, will be shown to be In our previous work we utilized this operator as the generator or the unitary component of the RE, for a ferromagnetic spin chain coupled by traceless operators to independent baths, with a time-independent system Hamiltonian. At that time we did not realize that it matched the Hamiltonian derived in Ref. [33]. Equation 16 generalizes the unitary generator to wider range of open quantum systems, within the RE framework.
To describe the dissipative dynamics of the RE, we introduce the dynamical tensor D(t), also in the interaction picture, as The relaxation tensor will be defined by the realignment of D(t) further below. The dynamical tensor is not necessary for calculations, since we can rewrite Eq. 14 in terms of (rank-2) matrices only, as after some algebra. The point of splitting the RE like this is that it disentangles the unitary and the dissipative dynamics, in a uniquely defined way. To fully grasp and then prove this statement we need to invoke the C algebra M (N ) of the N × N complex matrices next.

B. Superoperator Representation
Any trace preserving and Hermiticity preserving ME can be expressed in terms of the time-dependent generator of quantum dynamics [8], L t , which is a superoperator in M (N ), E ij = |i j| make the complete orthonormal set (COS) in M (N ). L(t) is Hermitian, [e.g., L ni,mj (t) = L mj,ni (t)] and satisfies n L ni,nj (t) = 0, ∀i, j, which assures trace preservation.
In the RE in particular, we assume for now that the system Hamiltonian is independent of time and |i is the energy eigenbasis. Thus we can apply the eigenoperator decomposition given by Eq. 8. Then, from Eq. 14 and after some algebra, we determine We can further recast the generator, as follows. First, using n L ni,nj (t) = 0, ∀i, j, as noted above, the RHS of Eq. 19 can be written down as (21) Second, let D(t) = ni,mj D nm,ij (t)E ni E † mj be the representation of the dynamical tensor (Eq. 17) in the eigenbasis. Define the relaxation tensor G(t) to be the realignment G ni,mj (t) = D nm,ij (t). G(t) is hermitian, e.g., G ni,mj (t) = G mj,ni (t), but D(t) is not generally so. Then, after some algebra Eq. 18 is recast as In third case, we can write the RE as where H ls (1) is given by the first line in Eq. 16, and Thus, in the basis {E ij }, there are at least three different ways to split the RE between the dissipativelike and the unitary-like terms, e.g., into terms like These forms are all GKSL-like, which demonstrates non-uniqueness of the dissipative-like and the unitarylike components of the RE. This implies that the dissipative-like terms generally carry information about the unitary dynamics, and vice versa.
Our physical idea of the relaxation tensor would be that it carries all the information there is about dephasing and other nonunitary processes in the system, but is not entangled with any unitary dynamics. To find the relaxation tensor, we need to identify which form, out of Eqs. 21-23, or any other, will not alter the unitary dynamics if we change the basis in M (N ). Such form is unique and has the Hamiltonian H ls (t) and the dynamical tensor D(t).
Proof: First we sketch the proof in Fig. 1 to gain conceptual understanding. Reference frame (a) represents M (N ) in the basis {E ij }, which does not include the unity matrix. Reference frame (b) also represents M (N ), but the basis now includes the unity matrix. Reference frames (a) and (b) are connected by a unitary transformation U . The generator of the RE can be split between the unitary-like and dissipative-like components, as we just discussed. The ellipses in (a) represent different splits of the generator, with dissipative-like and unitary-like component implied by circles and stars, re-spectively. In frame (b) the dissipative and unitary components are uniquely defined by the presence of the unity matrix in the basis [8,Lemma 2.3]. In frame (a), on the other hand, out of the many ways to split the Redfield generator, there can be only one with both the relaxation tensor and the Hamiltonian connected to the unique two in frame (b) via the unitary transformation. For this particular splitting, the dissipative and the unitary dynamics are disentangled since changing the basis in M (N ) changes neither the relaxation tensor nor the Hamiltonian. Thus, the goal is to transform the dissipative and the unitary generators from frame (b) to frame (a). Now we proceed to the mathematical execution.
We introduce a suitable basis in M (N ) and apply the proof of [8,Lemma 2.3]. Consider an operator in M (N ) with the properties trL(X) = 0 and [L(X)] † = LX † . As an example, L ni,mj could be the matrix L(t) from Eq. 19 at fixed time, but the following procedure can be applied to any trace-preserving and Hermiticity-preserving linear map. Now, switch the basis in M (N ) as follows: Importantly, F N N = (1/N ) 1/2 1, while the rest of the basis operators are traceless. Substituting, we find Here with primes indicating unequal indices. Now, [8, Lemma 2.3] applies to Eq. 27 because the basis includes the unity matrix, leading to the unique separation of the unitary and dissipative dynamics: The unique traceless Hamiltonian H is determined using Transforming back to the basis {E ij }, after some algebra, we obtain which corresponds to the traceless Lamb shift The remaining terms, e.g., L ni,mj (t), n = i, m = j; c ni,p (t), n = i, p = 1, ..., N − 1; c p,ni (t), n = i, p = 1, ..., N −1; and c pq (t), p, q = 1, ..., N −1 make the unique The interest in this matrix is that, for Markovian dynamics, its PSD is equivalent to the condition of CP of the map [8,9], while for non-Markovian dynamics, it is only a sufficient condition [36]. But here, we are primarily interested in its uniqueness property.
To determine the relaxation tensor in the frame {E ni }, we again change the basis from F ij to E ij . This is equivalent to applying the unitary transformation U † onto ma- while all other elements outside the block are zero.] Thus the relaxation tensor (viewed as an N 2 × N 2 matrix) and the Kossakowski matrix, respectively G and K, share the same eigenvalues, G has an extra eigenvalue of 0, and the PSD of one matrix is equivalent to that of the other. Since K is uniquely defined [8, Lemma 2.3], G will also be unique. After some algebra, the relaxation tensor can be shown to be (δ ni L kk,mj + L ni,kk δ jm ).
(34) The derived Hamiltonian and the relaxation tensor apply to any trace-preserving and Hermiticity preserving generator. Let us apply Eqs. 33 and 34 onto the coefficients given by Eq. 20. After tedious algebra and realignment, we obtain the traceless generator of the unitary dynamics and the dynamical tensor given by Eqs. 16 and 17, respectively. QED.
The extension to a time-dependent system Hamiltonian in the RE is straightforward but cumbersome, since the frequency integrals cannot be truncated by the Bohr frequencies. Nevertheless, Eqs. 16-18 are in the interaction picture and independent of the basis, and therefore apply to the time dependent system Hamiltonian as well.
Since the Kossakowski matrix and the relaxation tensor are essentially the same entity, we will use the terms interchangeably for the rest of the paper. (It will be implicitly understood that they refer to a different COS in M (N ), including the dimension difference.) If A α are all traceless, then Λ αβ (t) will also be traceless, and the Lamb shift and the dynamical tensor are given only by the first lines in Eqs. 16 and 17, respectively (the bottom lines in those equations do not affect the dynamics). We applied this simplified Lamb shift and the dynamical tensor in our prior work [14], based on a reasonable guess as to how to split the Redfield generator between the unitary and the dissipative component. Namely, we noticed that extracting this particular Hamiltonian out of the generator leaves the dissipative component amenable to further approximations. Independently from us, this particular splitting was also studied in Ref. [33].
In summary, within the second order perturbative RE, H ls (t) generates the unique unitary quantum dynamics in the interaction picture. The above proof can be utilized to obtain the Lamb shift and the dissipator at higher order in the coupling, in the time convolutionless master equation. The unitary component of the reduced dynamics is universal, since it is CP and well controlled by the perturbation theory, irrespective of what kind of environment the system is exposed to. Any additional approximation, for example to impose the CP condition on the dissipator, which we will study next, breaks this universality and should not be performed on the unitary generator. Since it will inherit the Lamb shift, such further approximated CP master equation will be less approximated relative to that where the unitary generator is approximated consistently with the dissipator.

II. THE GEOMETRIC-ARITHMETIC MASTER EQUATION
Now we proceed to impose the GKSL form on the dissipator of the RE, differently from how this was done in the singular limit or the SA. Now we allow for H S to have time dependence, in which case |i will be an arbitrary basis (otherwise, it will be the eigenbasis of H S ). For simplicity we assume that the coupling operators are all traceless, so only the top most line in Eq. 17 matters. Then we can rewrite the dynamical tensor as Since in the Markovian limit, the following result for the dissipator will be identical to that of the ULE, let us first examine the dynamics at that time, where Γ αβ (ω, t) is replaced by the asymptotic value leading to the dynamical tensor Next, we apply the crucial geometric-arithmetic mean approximation on the term Γ αβ (ω) + Γ βα (ω ), which is the essence of the GAME. Note that the frequencies of the system are expressed through the Fourier transform of the coupling operator, introducing the need of a frequency order ω c ω A , where ω c and ω A are the frequency ranges of the spectral densities of the bath and the system operators, respectively.
To derive the GAME, note that the relevant frequency range in the sum in Eq. 36 is given by |ω|, |ω | min(ω c , ω A ) = ω A . Just like in the SA, at time t the oscillations at frequency difference |ω−ω | 1/t will "average out". The essence of the matter is that at time t 1/ω A , the inequality |ω −ω | |ω|, |ω | will hold at the relevant frequencies. Thus, keeping in mind that J αβ (ω) is PSD, we can apply the matrix geometric-arithmetic mean approximation: At the relevant frequencies the error decreases inversely with time, that is, , so the relative error will be ∝ 1/ω A t. Thus, the region of applicability of the GAME, set by the domain of the geometricarithmetic mean approximation, is that the time of the dynamics be 1/ω A , which defines the corner time of the approximation. This is a much weaker condition than what is required for the validity of the SA, where the smallest oscillation frequency relegates the corner time. In other words, since ω A , being the typical frequency in A(t), is usually much higher than δ = min|ω| > 0, the GAME can be used even if the relaxation time t is 1/δ. The claim 2 from Refs. [15,31], that the ULE (or the GAME in our case) has a region of validity independent of the system frequency, is incorrect. Only if ω A < ω c such claim would be valid. Claim 2 was an unfortunate result of a rigorous but too loose of an error bound, which terminated any dependence of the error on the system frequencies. The dependence of the accuracy on the system frequencies will be explicitly studied and compared with the SA in section II A.
In Eq. 38, the square root is the matrix square root, and the ordering of frequencies is important. [Since J(ω) is PSD, so will be the square root; we could have done much better on the RHS, e.g., use [ J(ω ) J(ω) + J(ω) J(ω )] αβ instead, but unfortunately that will not lead to the GKSL form we seek.] Substituting Eq. 38 into Eq. 36 we find where Here we define the jump correlator as matrix which is the same as in the ULE. Now substitute Eq. 39 into Eq. 22, and the ME assumes the GKSL form, e.g., In this equation, M µ (t) represents M ∞µ (t), but it will be generalized shortly. For a time-independent Hamiltonian H S , the time dependence in M ∞µ (t) is entirely due to the interaction picture that we utilize, and in the Schrödinger picture the ME is time-independent consistent with Markovian dynamics. Next we consider the dynamics at time scale t τ c , where the GKSL form emerges as in Ref. [33]. Namely, Eq. 12 can be approximated at zeroth order in α as It is easily shown that matrix C αβ (0) is also PSD, and after similar algebra, we find where and the ME again assumes the GKSL form. At the very problematic intermediate time scale, defined in the case of the GAME as τ c < t < 1/ω A , we cannot rely on the geometric-arithmetic mean approximation. Short of anything better, we utilize the Jordan decomposition J(ω, t) = J + (ω, t) − J − (ω, t) and interpolate. Namely, since per the above the negative component J − (ω, t) approaches zero at t τ c and t 1/ω A , we interpolate by dropping J − (ω, t), so that Here we introduced the timed jump correlator to be the matrix The dependence on two time arguments stems from the memory in the bath. Fig. 2 displays the bath and the jump correlation functions, for the single Ohmic bath at zero temperature and spectral density where α is the dimensionless coupling constant, Θ(ω) is the Heaviside step function and ω c is the environmental cut-off frequency. The jump correlator has very complicated time dependence, and ranges from a delta at τ → 0 to that given by the long time limit in Eq. 52. The complicated time dependence reflects forcing of the GKSL form onto a master equation that naturally resists that form at the intermediate time scale. As Fig. 2(c) shows, the interpolation is akin to trying to approximate a sphere by a horse. We are curious if the recent regularization method from Ref. [30] can be adapted to better address this time range. For the Ohmic bath, the BCF is given by while the jump correlator in the long time limit from Eq. 49 is The timed spectral density according to Eq. 12 is explicitly A. GAME as an Advanced SA In the C algebra M (N ) a state ρ in the Schrödinger picture is replaced by a column vector |ρ made up of the columns of the matrix appended one after the other. In this format, the generator of the RE, L, is a superoperator with matrix elements L nm,ij = L ni,mj . The ME in the Schrödinger picture is then We assume that the master equation is asymptotically Markovian, e.g., L is time-independent. Since we operate within the weak coupling theory, the rate in the interaction picture (||˙ ||) is much smaller than the inverse bath correlation time.
On a time scale much shorter than 1/||˙ ||, the dynamics of the density matrix is approximately the same as the effect of coarse-graining in the interaction picture. To see this, note that the transformation of superoperators between the Schrödinger and the interaction picture is L int = e −L0τ Le L0τ , where L 0 = −i[H S , •] is the free system Liouvillian. The ME in the interaction picture is equivalent to the following: If t is much smaller than the inverse rate 1/||˙ || ∝ 1/α, then we will make the approximation (τ ) ≈ (0), so Eq. 55 becomes This equation shows the equivalence between the dynamics and the coarse-graining. The coarse-grained superoperator is defined as while the ergodic averagẽ is the superoperator in the SA. The SA becomes valid when t becomes higher than the highest oscillation period.
In practice, the coarse-graining amounts to replacing L ni,mj →L ni,mj = L ni,mj sinc[(ω nm − ω ij )t/2], where sinc(x) = sin(x)/x. (Note the useful property that both ω nm − ω ij and the Frobenius norm, that we use next, are invariant with respect to realignment ni, mj → nm, ij.) The difference between the solutions of the RE and the SA, as a function of time t, is tightly bound by the Frobenius distance ||L t,int −S||, while that between the RE and the GAME is bounded by ||L t,int −G t,int ||, whereG t,int is the coarse-grained superoperator of the GAME.
As an example, we consider a two-qubit system governed by the system Hamiltonian with the traceless coupling operator to an Ohmic environment The spectral properties of the environment are given by Eqs. 50, 51, and 53, and we study the dynamics at zero temperature. Here we assume the detuning δ = 0.01 and for the spectral density parameters we use ω c = 10∆.
The Frobenius distances decrease, on average, inversely with time, as shown in Fig. 3. The corner time for the error drop in the GAME, indicated by the red arrow, is the inverse drive 1/∆. This disconfirms the notion that the validity of the ULE is determined solely by the bath correlation time [15].
The blue line displays the Frobenius distance between the coarse-grained Redfield superoperator and the ergodic average. The corner time, indicated by the blue arrow, is approximately 1/δ, e.g., comparable to the longest oscillation time of the system. Thus, while both regions of validity are governed by the level structure of the system, the GAME has a much wider region than the SA. In particular, the GAME will be accurate if the relaxation time is much longer than 1/ω A , e.g., if Γ R ω A , where Γ R is the relaxation rate. In complex, many-body systems, that usually have exponentially small level spacings, this condition can be satisfied, in contrast to the SA where the requirement Γ R δ is effectively impossible.

III. MARKOVIAN DYNAMICS AS AN INTERTWINING MAP
In this section, we examine the asymptotic Markovian dynamics, outside the domain of the SA, and find that the exact dynamics will be governed by a GKSL equation if and only if the environment is in the singular coupling limit. In other words, there is nothing abnormal about the Markovian equation not preserving CP. Our proof is similar to that by Whitney [33], who determined that, in an Ohmic bath, the GKSL form of the RE can occur only in the extreme high temperature limit, when the environmental correlation time is negligibly small compared to all other relevant time scales.
In the singular coupling limit, C(t) ∝ δ(t), which is equivalent to the frequency independencies (flatness) of the spectral density. The "if" part is easy to prove: the RE in the singular coupling limit has the GKSL form, the fourth order cumulant is zero [46, Eqs. B3 and B4], and more generally, the exact ME is also in the GKLS form [47].
We proceed to the "only if" part. First we demonstrate that the asymptotic RE violates the GKSL form in proportion with the variance in the spectral density. For simplicity we assume single heat bath as in Ref. [33], so that we can drop indices α and β in Eqs. 16 and 17. Also, here we assume that the system Hamiltonian is time-independent, and take advantage of the eigenoperator decomposition 8. Using the identity n|A(ω)|i = [A(ω)] ni = A ni δ ω,ωin , where ω in = E i − E n are the Bohr frequencies, we find the relaxation tensor for traceless coupling operator, by realigning Eq. 17, We introduce single indices in lexicographic order, p = in, make the time dependence implicit, and simplify the notation as Γ(ω in , t) = Γ p = J p + iS p . Then we solve the eigenvalue problem for G (see Appendix B). As in Ref. [33], there are only two nonzero eigenvalues, where J = p |A p | 2 J p / p |A p | 2 , Var(J) and Var(S) are the respective variances, mutadis mutandis, and || || is the Frobenius norm. The eigenvalues are proportional to α and λ − < 0. One exception is that of λ − = 0, iff Var(J) = Var(S) = 0, i.e., in the singular coupling limit. Next examine these eigenvalues in the case of the Ohmic bath, with the BCF 51. Fig. 4(a) displays the ratio of λ − /λ + , for the system in Sec. II A at δ = (1/2)∆, and for the RE with time-dependent coefficients. At t → ∞, λ − /λ + ≈ −1/3, which could be misinterpreted as significant nonpositivity, and hence significant CP violation.
In Fig. 4(a) we also plot the ratio of the sum of the negative and the positive eigenvalues of the relaxation tensor in the TCL4 t ME, and observe that the nonpositivity is not cancelled by the increasing perturbative accuracy. In fact, the ratio becomes slightly more negative [see the inset in Fig. 4(a)]. This disconfirms the proposal in Ref. [32] that involving the fourth order terms can fix this nonpositivity problem, since the Markovian ME naturally does not preserve CP and needs no fixing, if it works adequately. Now we generalize this analysis to the exact TCL equation, where the relaxation tensor G can be expressed perturbatively as and k p=1 α p G (2p) (t) is that tensor in the TCL2k t ME. By the perturbation theory, the eigenvalues of G are Since λ − is proportional to α, the eigenvalue λ − + O(α 2 ) remains negative for sufficiently small α. Thus, the exact asymptotic TCL ME is not in the GKSL form. Now, the question is how will the reduced dynamics be CP, if the exact relaxation tensor is perpetually nonpositive? In the pre-Markovian stage where the ME is time-dependent, the answer is that the PSD of the timedependent Kossakowski matrix is sufficient, but not necessary condition for the CP of the map; see Ref. [36] for examples of this and the discussion of the related intertwining maps. The latter are composite maps: the states of the system at time t 1 > 0 are sent to the initial factorized state, by the inverse map (assuming that it exists), which is then sent to the state at time t 2 > t 1 by the time-forward map. Since the Markovian dynamics is asymptotic, it is an intertwining map in the limit t 1 1/ω c and t 2 > t 1 . The intertwining map need not be CP, while the quantum map initiated into the factorized state is of course CP.
To demonstrate that this is indeed the case, we invoke the TCL4 t master equation, and show that it exhibits a strong improvement toward CP of the quantum dynamical map, in absence of any suppression of the negative eigenvalues of the relaxation tensor, confirming that the non-CP Markovian dynamics intertwines a completely positive map. The measure of CP is obtained by utilizing the Choi matrix [48], defined here as It is easily verified that tr(C t ) = 1. The dynamical map generated by L t will be CP, iff C t is PSD [48].
In the example in Fig. 4, we numerically solve the ME for every matrix E ij as the initial condition, and then compute the Choi matrix. The violation of CP of the map is measured here by the sum of the negative eigenvalues of the matrix (Σ − ) 1 We note in passing, that if the steady in the Schrödinger picture [ρ(∞)] is unique, one can show So, the PSD of the asymptotic state is equivalent to the CP of the injection of all density matrices in M (N ) into ρ(∞). Fig. 4(b) displays Σ − versus time. In contrast to the weak dependence of λ − /λ + on the order in α [ Fig. 4(a)], increasing that order strongly suppresses the nonpositivity of the dynamical map.
In Fig. 4(c), we display the non-PSD of the Choi matrix as a function of α. Σ − is the extremal value between t = 0 and t → ∞. For n = 1 and 2, it scales as α n , where 1 As a calibration, if there is no CP violation, then Σ − = 0 and the sum of the positive eigenvalues (Σ + ) will be one. If there is a CP violation, Σ − < 0. Since Σ − + Σ + = tr(Ct) = 1, the relative change in Σ + will not be significant, if |Σ − | 1. Thus, the characteristic scale of |Σ − | is one: the smaller the |Σ − | relative to one, the weaker the CP violation. 2n is the perturbation order. As n increases, quantum dynamics approaches CP while the master equation does not approach the GKSL form. Thus, the non-CP Markovian dynamics intertwines a non-Markovian CP map. There is a another very practical significance to determining the Choi matrix. The strong suppression of |Σ − | between the RE t and TCL4 t cannot be by accident. TCL2n t master equations are very cumbersome, making it easy to make a mistake. Thus, if we confirm that Σ − scales as α n , this will imply that all the nonpositivities at lower orders in α are properly cancelled out, which will assure the correctness of the calculation of the TCL2n t . Thus, we can confirm the correctness without knowing the exact solution.

IV. SYSTEMATIC ACCURACY OF THE GAME
Here we expose and evaluate the non-universality of the GAME numerically. On the contrary, the RE is universal in a sense that no matter what the state of the system is, both the coherence norms and the population norms, as defined in the introduction, will be accurate in the leading order of α throughout the full course of the dynamics. In the steady state of the RE, for instance, it is known that the populations and the coherences are of zeroth and first order in α, respectively [32,[49][50][51]. By ) ) ) ) ) ) ) ) ) . & ! % ) ) ) ) ) ) ) ) ) ) the universality, the respective errors will be of first and second order.
Here we will find that if the initial state has coherences at first order in α, then during the entire dynamics the coherences will remain at first order in α. This dynamics, where the coherences are at first order in α, is what we mean by incoherent dynamics. As an example, the dynamics of the system will be incoherent if initiated to an excited energy eigenstate. The RE is universal in a sense that the errors of the coherences are at second order in α throughout this dynamics.
Any attempt to enforce complete positivity onto the incoherent dynamics will result in loss of accuracy in the coherences, downgrading the error to be at first order in α. That means that the relative error of the coherences will be at zeroth order in α, which is non-universal. This breakdown of universality has been studied in the steady states in Ref. [32]. Here, we extend this result to all incoherent dynamics, by examples.
Remarkably, coherent quantum dynamics, where the coherences are at zeroth order in α, is systematically accurate when modeled by the GAME. This is the main result of this paper. Thus, if we ask that the GAME describes the density matrix with both the populations and the coherences at zeroth order in α, then both errors will be at first order in α. Although not true in general, in many interesting physical situations such condition is valid.
First we study the example from the previous section, in strong detuning: δ = ∆/2. Initially we examine the situation where the qubits are prepared in the highest energy eigenstate, |1, 1 , insuring incoherent dynamics. We determine the ME solutions using the Runge Kutta 4 method. The relaxation rate (Γ R ) is determined by fitting the population of the initial state versus time to the exponential, leading to Γ R = 1.084α∆.
) ) ) ) ) ) ) ) ) . ( ' ) ) ) ) ) ) ) ) ) Initially we investigate the accuracies of the MEs in the Markovian limit, by comparing the states obtained by solving the ULE and the GAME with those from the RE. Here we define the state error by the trace distance to the RE state. Crucially, we separate the diagonal and offdiagonal parts, i.e., ρ = ρ d + ρ c , and measure the respective errors in the population and coherence separately, as defined in the introduction. This is different from previous practice where the accuracies were quantified by the trace distance between the full states [14,15,25]. The reason to separate the errors is that the RE is known to have reduced absolute (but not relative) accuracy of the populations [32,[49][50][51]. 2 The top green line in Fig. 5 displays the coherence of the state versus time obtained from the RE. The average coherence in this time interval is 3.26α. The dottedblue and the dashed-red thin lines display the errors of the coherences of the ULE and the GAME, respectively. The error in the GAME is 0.63 times that of the ULE, confirming the supremacy of the GAME, e.g., that the universal Lamb shift leads to more accurate results. In both equations, the relative error is at zeroth order in α, demonstrating their non-universality. The average error in the coherence of the GAME is 0.412α, or approximately 12.6%. Although not universal, in many interesting situations the error is good for qualitative description. On the other hand, the errors of the populations in the GAME and the ULE are systematic, e.g., at first order in α, and virtually the same. Now we discuss the accuracies of the non-Markovian master equations. Namely, we evaluate the trace distances between the RE t and the GAME t with respect to : the TCL4 t solution, which has negligible error in comparison. The results are displayed in Fig. 6. The coherence of the TCL4 t (green line) has the average of 2.95α and the trace distance between the coherence of the GAME t and the TCL4 t (red-dashed line) are very similar to the corresponding curves in Fig. 5. Similar can be said about the accuracies of the populations. There is, however, a dramatic enhancement in the accuracy of the coherence of the RE t equation, as shown by the black-dotted line in Fig. 6. Namely, the error in the coherence of the RE t is factor of O(α) smaller relative to that obtained by the GAME t . Fig. 7 displays the average coherence and the error as a function of α. The average coherence is obtained by averaging on the time interval [0, τ R ], where τ R = 1/Γ R . In Fig. 7(a), the average coherence (determined from the TCL4 t ME) scales proportionally with α, as shown by the green circles. The accuracy of the coherence determined by the GAME t (red triangles) also scales proportionally with α. Thus, the relative accuracy is independent of α, implying that the GAME t is not universal. By contrast, the accuracy of the coherence obtained by the RE t scales quadratically with α. Now, as elaborated in Ref. [32], if the figure of merit is the trace distance between the full states, then the error will be governed by the population and one can mistakenly claim that the RE t and the GAME t are on equivalent level of accuracy, as was done in the ULE paper [15]. The upshot is that one must compare coherences to coherences and populations to populations, which exposes the non-universality of the GKSL equation.
Next we discuss the numerical results on the relaxation from an initially superimposed system state (|01 − |10 )/ √ 2, and show the results for the coherence in Fig. 7(b). The average coherence in the time interval [1, τ R ] is comparable to the initial coherence of 0.5, since the dephasing time is long (e.g., T 2 = 2τ R ). Against the backdrop of strong coherence, the error in the coherence of both the RE t and the GAME t scale proportionally with α. In this sense the systematic error of the GAME is restored for coherent quantum dynamics.

A. In-Depth Study: Heisenberg Spin-Chain
Next we examine the accuracy of the GAME in a significantly more complex system and evaluate how the systematic accuracy of the GAME breaks down as the dynamics approaches the incoherent stage. We consider a spin-1/2 ferromagnetic chain with N = 19 sites, interacting with nearest neighbour exchange interaction and uniaxial anisotropy due to dipole-dipole coupling: We introduced this model in our prior work [14], to study the accuracy of the GAME. Note that [H S , S z ] = 0, so that S z is a good quantum number. The characteristic frequency of the system is given by the energy gap ∆, which corresponds to the ferromagnetic resonance frequency. The eigenenergies of H S versus S z are displayed in Fig. 8(a). The model exhibits a wide distribution of frequencies associated with spin modes, including ferromagnetic resonance ∆, standing spin-waves, and domain walls [14]. The spin chain is coupled to 3N independent Ohmic baths with the BCF given by Eq. 51 and system coupling operators s x , s y and s z at each site. The range of relaxation rates is very broad, from that of the ferromagnetic resonance (∝ α∆) to the very slow motion of domain walls near the center of the chain [14].
In the previous work, we determined the error of the GAME by taking the trace distance between the solutions of the RE and the GAME. But here, we compare coherences to coherences, since this is where the systematic accuracy can fail, as discussed in the previous example. We have checked that the populations always have errors at first order in α, thus we will not discuss the populations. We proceed to map out the errors of the coherences as a function of the state, to determine the range of states where the GAME is systematically accurate.
The system is prepared into the initial state with magnetization along x axis, perpendicular to the chain. That is, the initial state is the eigenstate with S x = 9.5. To study incoherent dynamics, we will set all the initial coherences to zero, while the populations will be the same as in this eigenstate. Fig. 9 displays the dependence of the coherence of the RE states and the error of the GAME coherences, e.g., (1/2)||ρ RE,c || 1 and (1/2)||ρ GAM E,c − ρ RE,c || 1 , respectively, versus scaled time t/T 1 , where T 1 is defined in Fig. 8(b). Fig. 9(a) displays coherent dynamics with ρ(0) = |S x = 9/2 S x = 9/2|. The top three lines display the coherences for α = 0.0001, 0.001 and 0.01. The physical time range is the same between the curves, which is why the ranges differ in the scaled time. We note that the coherences start to drop at t/T 1 < 1, because T 2 is smaller than T 1 [ Fig. 8(b)].
The errors of the coherences are displayed by the corresponding bottom three lines, indicated by the boxed values of α. Initially, (1/2)||ρ GAM E,c −ρ RE,c || 1 increases linearly with time but saturates at the time scale of ∝ 1/∆, at the corner times indicated by the crosses in Fig. 9(a). The crosses occur at the same physical time, but are displaced in the scaled time t/T 1 . These results confirm the same behavior as we studied in the trace distance of the full states in the previous work [14]. At saturation the coherence errors are proportional to α, e.g., roughly 10α.
The key observation is in the following: As the coherences begin to drop versus increasing time, so do the errors. In other words, the relative error of the coherence is independent of the coherence. In this example, ||ρ GAM E,c − ρ RE,c || 1 ≈ 10α||ρ RE,c || 1 . Thus, the error is systematic.
The systematic regime persists until the coherences drop below a value indicated by the stars in Fig. 9(a) (the red star is only suggestive). For α = 0.01, with further drop in the coherence the relative error also drops until it becomes comparable to α 0 . Thus, we conclude that the systematic regime will be valid if the coherences are above a threshold proportional to α.
We present incoherent dynamics shown in Fig. 9(b), where ρ(0) = diag(|S x = 9/2 S x = 9/2|). For each color, the top curve is the coherence and the bottom curve is the error of the coherence. These curves saturate the same way as the errors in Fig. 9(a). We observe in (b) that the coherences and the errors of the coherences are at the same order in α. So this is nonsystematic regime of the GAME. If we compare the coherent and incoher- The conclusion is that the GAME has relative error of the coherences at first order in α. But this fails when the coherences are on the order of α. Since α is usually very small in quantum technologies, this shows that the GAME can be used with systematic accuracy for a large range of quantum states and quantum systems. This reconciles the opposing claims of universality of the GKSL ULE [15] and the fallacies of the GKSL enforced RE [32].

V. EXAMPLE: DYNAMICAL DECOUPLING
As we have seen, the GKSL form of a master equation is natural only in the singular limit of the BCF and in the SA. For master equations not in that form, schemes have been developed to find an appropriate GKSL-form. Due to the lack of perturbative rigor in those schemes, it becomes a natural question why the enforcement should be done if at all. Here we show that enforcing the GKSL form is the only way to treat the problem of dynamical decoupling (DD) in the ME framework.
The DD protocol involves applying a time-dependent system Hamiltonian on a quantum system, in the form of a sequence of pulses [52,53]. In Markovian dynamics, the DD has been studied in Refs. [25,54,55]. If the time between the pulses is shorter than the environmental correlation time, it can cause dephasing suppression which may be important in quantum error correction [56,57].
Here we study an example of a spin-1/2 boson problem in the regime of pure dephasing. For simplicity we assume that the free Hamiltonian is zero, so that the system Hamiltonian is solely due to the external field, The spin-1/2 is flipped by a π pulse periodically with a period τ , and is coupled to Ohmic bath at zero temperature with spectral density given by Eq. 50, with traceless operator A = σ x /2. In the interaction picture, A(t) undergoes a series of discrete sign changes at times t = nτ . An alternative view of this, based on Eq. 15, is that A(t) does not have discrete sign changes, (but keeps the continuous unitary changes, if any), while the BCF absorbs the sign changes, as shown by the dashed red line in Fig. 2(a). Then, the effect of the pulses will be twofold. First, the effective spectral density, given by the Fourier transform of the apparent BCF, will be strongly suppressed if τ 1/ω c . Second, the correlation time of the bath becomes effectively τ c , making it more singular-coupling like.
First we solve the DD problem by the path integral technique known as Time Evolving Matrix Product Operator (TEMPO). The introduction to the method and the python code we use are available in Ref. [58]. The TEMPO algorithm can also be modified to calculate the process tensor [59]. This alternative formulation has been used in both optimizing quantum control procedures [60] and in calculating exact bath dynamics [61]. An open source python package [62] improves the performance and includes new approaches. Here we work in the Schrödinger picture. In the DD, the unitary propagator imposed by a single pulse will be U 0 ⊗ 1, where U 0 = −iσ y . Let the unitary propagator for the complete system, in time interval τ between the pulses, be W τ . Then, since the unitary dynamics is divisible, the unitary propagator at time t = N τ •] is the Liouvillian for a single pulse, while L τ = e −i[H B +H I ,•]τ is the Liouvillian for the continuous evolution of the system and bath between consecutive pulses.
Then, the reduced state of the system after N pulses in the Schrödinger picture is Next, the state can be written as a path sum over system states, by inserting identity resolutions between each L 0 and L τ . Due to the form of the Liouvillian L τ , tracing over the bath can be accounted for by a discretized Feynman-Vernon functional with N time steps [63]. An important simplification here is that there is no need to subdivide the time step of length τ and use a Trotter splitting, because the system Hamiltonian (equal to zero) commutes with H B + H I between the pulses. TEMPO utilizes matrix product states to store the Nindex tensor of the system, known as the augmented density tensor, which is the central entity in the exact Quasi Adiabatic Path Integral method (QUAPI) [64][65][66][67]. In TEMPO, the tensor is efficiently compressed by singular value decomposition (SVD), thereby making the calculation of the path sum numerically tractable. Out of the three errors in TEMPO, namely, the Trotter splitting error, the finite memory error, and the SVD truncation error, our application only exhibits the last one. We calculate using the open source code available in Ref. [58], except that we made a minor modification to take advantage of the absence of the Trotter splitting.
The degenerate spin-1/2 is initialized into |1 with spin along +z axis, and interacts with the bath and the ex-ternal pulses. In the asymptotic Markovian limit, the precise position of the pulse relative to the initial time has a weak effect and we shall not discuss it here. The result for the population as a function of time is shown in Fig. 10(a) and (c), for different pulse lengths. Each curve has two branches reflecting the fact that s z is flipped half of the time. Decreasing the pulse length below the correlation time 1/ω c causes a rapid increase in the dephasing time T 2 .
Next we solve the same problem using the GAME t master equation and show the results in Fig. 10(b) and (d), indicating remarkable agreement with the exact solution. By fitting a branch to an exponential, we obtain the dephasing time versus τ pulse , and compare the results from TEMPO and the GAME t in Fig. 11(a). The agreement is great. This is not surprising, because the DD protocol suppresses the effective correlation time in the environment. Now, if we use the RE t , the population versus time becomes nonsensical, as shown in in Fig. 11(b). This is due to the instability of the RE in the absence of the positivity guarantee. Thus, in this example the enforcement of the GKSL form is not only highly accurate, but also necessary.

VI. DISCUSSION AND CONCLUSION
In this paper we studied the nonpreservation of CP in perturbative, quantum master equations. Due to the finite correlation time of the bath correlation function, Markovian dynamics is an intertwining map of a non-Markovian dynamics, and need not be CP as such. We explicitly demonstrate this by determining the Choi matrix of the non-Markovian master equation, which confirms the (near) CP-preservation, when the Markovian master equation is not CP-preserving.
We rederive the completely positive geometric arithmetic master equation in the spirit of the secular approximation. Rather than taking the ergodic average, we average the dynamics on time scale that permits replacing the arithmetic mean at various frequencies with the geometric mean, thereby transforming the Redfield equation into a GKSL master equation without sacrificing the nonsecular terms. The characteristic time scale that determines the accuracy of the GAME is given by the range of the spectral density of the system-bath coupling operator, rather than the minimum oscillation frequency of the system. The result is that the GAME can be used if the relaxation time is much smaller than the longest oscillation time of the system.
As a side result, we extracted the usefull Lamb shift and the Kossakowski matrix for an arbitrary Hermiticity preserving and trace preserving master equation, by changing the basis in the algebra of complex matrices.
We investigated the range of system states for which the GAME has systematic accuracy in the coupling rate to the bath. The coherent quantum dynamics can be calculated at the leading order in the coupling rate, but the incoherent dynamics cannot.
In the regime of dynamical decoupling, we find that the GAME becomes nearly exact, while the perturbative master equation becomes unstable. This reversal of fortune is due to the fact that the dynamical decoupling decreases the apparent bath correlation time, making the environment effectively singular, thereby enforcing the GKSL form on the master equation.
Further work will study the bath states and seek to quantify the entanglement between the system and the bath, to learn how and when it forms, and how this is all quantitatively related to the negative eigenvalues of the system's Kossakowski matrix. Path integral formulation TEMPO can determine the exact bath dynamics [61]. Since the TCL and the path integrals are equivalent [68], the irrelevant part of the density matrix in the former can potentially be determined using this technique, suggesting that the full density matrix can be determined. Alternative techniques to get insight into the bath could involve reaction coordinates, [69,70] auxiliary modes [71][72][73][74], counting statistics [75][76][77] and density matrix renormalization group [78][79][80][81].
The author would like to thank Donghyun Jin, Grihith Manchanda, Huo Chen, Nicholas Ezzell, Evgeny Mozgunov, Daniel Lidar, Archak Purkayastha, Gerald Fux and Brendon Lovett for discussion; and Gerald Fux for adding a feature to the open source package on our request. This research was supported by DOE contract DE-FG02-06ER46281. Additional support from the Georgia Tech Quantum Alliance (GTQA), a center funded by the Georgia Tech Institute of Electronics and Nanotechnology was used to develop exact numerical methods on entangled low-dimensional magnetic system. Here the short-hand notion meanŝ i = A(t i ), i . = 0, 1, 2, 3, t 0 = t ij = C(t i − t j ), i, j . = 0, 1, 2, 3.
Here we define four bath correlation functions where Γ jm (t) = Γ(ω jm , t) is the timed spectral density (See Eq. 53). Since the four correlation functions embody four point correlations in time, there is no escape from rank-4 tensors, in contrast to the RE which can be studied entirely in terms of (rank-2) matrices. Note that if the spectral density is flat or C(t) = cδ(t), then In sections III and IV, we determine the integrals using the Simpson rule. The smallest sub-step in time is 0.0025/∆, and increases by factor of two in each consecutive integral. Fig. 12 displays typical time dependence of a matrix element in D(t). The inset shows that the approach to the asymptotic Markovian dynamics is very slow, making it very time-consuming to determine the asymptotic limit. To extrapolate to infinity, we determine D ∞ as the average of D(t) over the last two or three periods. The period is determined as the time between successive major maxima in the matrix element versus time, and is found to be independent of the matrix ele-ment. The relative error in the average is estimated to be in the 10 −8 range. The Frobenius distance between D(t) and D ∞ decreases as inverse time, as shown in Fig. 12(b), suggesting that the Markovian limit exists, indeed.