Completely Positive, Simple, and Possibly Highly Accurate Approximation of the Redfield Equation

Here we present a Markovian master equation that approximates the Redfield equation, a well known non-Markovian master equation derived from first principles, without significantly compromising the range of applicability of the Redfield equation. Instead of full-scale coarse-graining, this approximation only truncates terms in the Redfield equation that average out over a time-scale typical of the heat bath. The key to this approximation is to properly renormalize the system Hamiltonian, followed by the Markovian arithmetic geometric mean approximation (MAGMA), that restores complete positivity of the master equation. Three applications of the Markovian master equation are presented and demonstrate its simplicity and accuracy. In an exactly solvable quantum dynamical problem of a three-level system, we find that the error of the approximate state is almost an order of magnitude lower than that obtained by solving the coarse-grained stochastic master equation. In the example of a four-level system, we find that the Markovian master equation is more accurate than the non-Markovian Redfield equation that it emulates. It is also shown that the Markovian master equation retains a near Born-Markov accuracy when applied to a complex many-body system, a ferromagnetic spin-1/2 chain with long-range, dipole-dipole interactions between 25 spins.


I. Introduction
Quantum correlations between particles in quantum systems are counterintuitive, giving us exotic algorithms that can greatly accelerate some calculations using quantum computing as opposed to conventional computing. It remains an unsettled issue, however, to describe correlated quantum dynamics in large many-body quantum systems with the precision required for quantum computing if the systems are coupled to an environment. In the Markov approximation, the time dependence of the reduced state of the system, can be approximated by a first order linear differential equation, known as the master equation, which greatly simplifies this description. 1 There is only one form of the master equation that is a completely positive map, the deeply respected Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) form of master equation. [2][3][4] The GKLS equation is a consequence of the Diracvon Neumann axioms, while the physical scope of the equation rests in the coarse-graining of the reduced system dynamics over a timescale much longer than the environment correlation time, τ c . 1 The time in which the environment changes the system state in the interaction picture, τ r , must therefore be much longer than τ c , so that the master equation can be useful. The condition τ c ≪ τ r is the Born-Markov approximation. 5 It has been very difficult to derive the GKLS equation from first principle calculations using the total Hamiltonian of the system and environment as a starting point. 6 The operator-projection technique has been the standard way to derive it, 4,[7][8][9] and was first done by Redfield. 10 Although the Redfield master equation is derived by applying the Born-Markov approximation, the equation is non-Markovian with a flaw: it does not guarantee positivity of the state. 6 This can lead to negative probabilities of observable, which violates the mathematical definition of probability. Negative states also pose problems in terms of how to quantify the entanglement of the system, which is important for quantum computing. One measure of entanglement is the negativity of the state's partial transpose, 11 which raises the question: Can we trust negativity of the partial transpose as genuine measure of entanglement if the approximate state itself is negative?
Although it has this problem of producing negative states, the Redfield equation has major advantages as well. The Redfield equation is a highly accurate master equation that is simple to apply and not axiomatic. Thus, it has found significant use in quantum chemistry. 12 Nevertheless, the issue of negative states remains unsettling and debated.
There have been several attempts to cure the issue of negative states by modifying the first prin-ciple equation into the GKLS form. This was done for the first time by Davies,13 who gave us the rotating-wave approximation (RWA), also known as the secular approximation. The approximation amounts to coarse-graining of the reduced state on a timescale comparable to the Heisenberg time. Since this time scales exponentially with system size, the RWA does not work well for large many-body quantum systems that generally have exponentially suppressed, level spacings with increasing size of the system. More recently, coarse-grained Markovian master equations have been derived from first principles [14][15][16] capable of capturing correlation effects that the RWA cannot. As an alternative to coarse-graining, reasonable guesses of completely positive master equations have also worked well. 17,18 In this paper, we obtain a completely positive quantum master equation by applying the Markovian arithmetic geometric mean approximation (MAGMA) on the spectral density after the appropriate renormalization of the system Hamiltonian. We identify and drop some terms in the Redfield equation that average out over typical timescales of the bath. We label these terms "detuning terms." Without coarse-graining the entire equation, we arrive at the Markovian equation.
We estimate the approximate state error bound and show that it is comparable to that of the coarse-grained stochastic equation (CGSE). 14  the exact quantum dynamics. We also explain how the error gets reduced. Lastly, we investigate a large many body interacting system, a ferromagnetic spin chain, and find that the equation remains accurate within an error close to the Born-Markov error.
The paper is organized as follows. In Sec. II we introduce various master equations relevant for this work. We also introduce the notation used similar to that of the recent paper by Mozgunov and Lidar, 16 but slightly compressed, namely, we use the Hadamard matrix product instead of the operators A ω in Refs. 4,14,16 . In Sec. II A, II B, and II C we review the Redfield, Coarse-Grained Redfield, and the Davies-Lindblad master equations, respectively, as a prelude to MAGMA, which is presented in Sec. II D. In Sec. III we implement MAGMA on an exactly solvable model, examine the error of the approximation, and give a new way of looking into correlated quantum dynamics in the presence of the heat bath. In Sec. IV we introduce a 4-level system that displays negative Redfield dynamics and examine the source of negative eigenvalues of the state and the physical aspect of restoring complete positivity. We also examine how restoring positivity decreases the error relative to the exact quantum state. We end the paper with an investigation of a 25-body quantum system in Sec. IV.

II. Introducing the problem
Here we consider a quantum system coupled to a heat-bath, represented by the total Hamiltonian H 0 , I 0 , and A are the respective Hamiltonian, unity, and arbitrary operators performing in the Hilbert space of the quantum system, H B , I B , and B are the corresponding operators for the heat bath, and A ⊗ B is the interaction Hamiltonian between the system and the heat bath. All the operators are assumed to be Hermitian. For notational clarity we only consider one coupling term in the interaction Hamiltonian. The results can be straightforwardly extended to a sum of interaction terms assuming uncorrelated heat baths.
In the Hilbert space of the quantum system, we use the eigenbasis of H 0 , H 0 |n = E n |n , where |n are the eigenvectors and E n are the eigenenergies. Then, in the interaction picture where A nm = n|A|m , ω mn = E m − E n , t is time, and we useh = 1. Operator A(t) can be represented in the matrix form as where • is the Hadamard product (element-wise product) and [Ω(t)] nm = e −iωmnt . Here we assume no degeneracies, so that, We introduce the bath correlation function and ρ B is the reduced density matrix of the heat bath. We assume a stationary heat bath so that [ρ B , H B ] = 0, and in that case C(t) = C ⋆ (−t). If necessary, see Ref. 16 for more details, but the information provided should be sufficient to follow the paper. Here we reserve the symbols ρ and ̺, for the density matrix in the Schrödinger and interaction picture, respectively.
, and g is the dimensionless system-bath coupling constant.

A. Redfield equation
The equation is a widely used fundamental master equation describing fluctuations and dissipation in quantal systems. 9 It is derived from first principles, by applying the Born-Markov approximation and projection operator technique. We begin here from Eq. 16 from the recent where is termed the "filtered" operator. 16 After calculating the integral, one finds A f = A • F , where F = 1 2 γ − iS, and matrices γ and S have elements γ nm = γ(ω nm ) and S nm = S(ω nm ), respectively. Function γ(x) is the spectral density, equal to the Fourier transform of the bath correlation function, and is non-negative. S(x) is equal to the Cauchy principal value on the half-range Fourier transform of the correlation function, and here we refer to it as the principal density. Table I summarizes key formulae for the bath correlation functions and its integral transforms. A specific example of the Ohmic heat-bath at zero temperature is included in the bottom line.
The range of applicability of the Redfield equation is given by the condition of the validity of the Born-Markov approximation. The condition can be obtained by estimating the fourth order term in the Dyson expansion on the system state, and comparing it to the second order term. 19,20 The approximation fails when the two terms are comparable. That leads to the range of validity given by τ c ≪ τ r , where τ r is the relaxation time of the reduced state in the interaction picture, and τ c is the timescale on which the correlations in the heat bath decay.
Rigorous error bound for the Born-Markov approximation exists. 13 Recently, it has been shown that the solution of the Redfield equation has the error bound of where τ sb is the lower bound of τ r (also known as the smallest relaxation time). 16 This error bound is general, and applies to any total Hamiltonian and any density matrix.
Here 1/2||ρ B (t) − ρ exact (t)|| 1 is the trace distance. (For a matrix X, the trace norm is defined as ||X|| 1 = T r √ X † X which is equal to the sum of the singular values.) The O symbol means "on the order off", and its precise definition is given in Ref. 16 .
The actual dynamics of the system depends on the total Hamiltonian and the initial density matrix, and the time scale τ r is usually larger than the bound τ sb . The state error should similarly be smaller than the bound given in Eq. 6. It is sensible to write the error as this, which is consistent with the error bound Eq. 6, (since τ r ≥ τ sb ), and approaches 1 when τ c ≈ τ r , as it should, because that is where the Born-Markov approximation fails. 19,20 Eq. 7 can be recast in a different form, With the renormalization of the system Hamiltonian, the Redfield equation becomes or in terms of matrix elements, Note that all the summands involve this one and only function of system frequencies, In contrast to the initial Eq. 4, in this form, the equation exhibits balanced state gains and losses, due to the coupling to the environment. This will be critical to our approximation, because it allows us to apply the same approximation uniformly on each term of the master equation.
The matrix elements of the renormalized Hamiltonian are The renormalization involves both principal and real densities S and γ, respectively, in contrast to the Lamb-shift in the RWA, e.g., see Eq. 24.
The Redfield equation in form 10 has been studied recently in context of embedding of the quantum system into an expanded system that displays Markovian quantum dynamics. 21 However, the density matrix of the subsystem, is presented by the off-diagonal block of the expanded state, and such embedding does not resolve negativity of the reduced state.

B. Coarse-grained Redfield Equation
Here we coarse-grain Eq. 11, as follows. First, we transform to the interaction picture [̺ nm = exp(iω nm t)ρ nm ], where the Redfield equation reads Next, we change the time symbol from t to τ , and time average on a Gaussian window. That is, we apply the convolution operator to the equation, while keeping ̺ at time t. T 0 is the coarse-graining time. This leads to whereG andH In the third step, we transform out of the interaction picture, [ρ nm = exp(−iω nm t)̺ nm ], and obtain the coarse-grained Redfield equation, The matrix elements of the renormalized Hamiltonian are Coarse-graining increases the error of the approximate state. Rigorous error bound exists for the state error added by the coarse-graining on a rectangular window, 16 where, we repeat, τ sb is the lower bound of τ r , and the LHS is the added coarse-graining error.
We use the Gaussian window instead of the rectangular one because the former makes it much easier to find this error. We obtain the same result as Eq. 21 but do not give proof since the proof essentially repeats that on the rectangular window. 16 Here it is sensible to assume that the actual error added by the coarse-graining is lower than the error bound given by Eq. 21 because the state evolves on the timescale τ r ≥ τ sb . In other words, to get the rigorous error bound, the coarse-graining time is compared to the minimum value of τ r (Eq. 21). On the other hand, to get the estimate of the error itself, the time should be compared to τ r , e.g., In other words, the replacement of τ sb in Eq. 21 with τ r , converts the error bound into the error estimate. This follows the same logic that we use to arrive at Eq. 7.

C. Rotating wave and CGSE approximations
This approximation amounts to discarding all but the following terms in Eq 11: Here, ǫ 0 nm = ǫ 0 n − ǫ 0 m , where ǫ 0 n = E n + δE L,n and is known as the Lamb shift.
The RWA can be re-derived directly from Eqs. 17, 19 and 20 by assuming T 0 ≫ 1/|ω nm |. In that limit, the exponentials can only be zero or one. As a result, the RWA is only valid when the coarse-graining time is longer than the Heisenberg time τ h = 1/δE, where δE is the smallest level spacing in the system. In that respect, the system relaxation time also needs to be of that order or longer, so that the coarse-grained equation can be of value. Therefore, the main disadvantage to the RWA is in its highly limited range of applicability. Note that at zero temperature, the dissipative part of the system dynamics under the RWA only involves positive frequencies (since the spectral density at negative frequency is zero). 22 More recently, a coarse-grained stochastic equation was derived 14 with better error and range of applicability than that of Davies', but not as low as the Redfield error bound given by Eq. 6.

D. A Markovian arithmetic geometric mean approximation (MAGMA)
Since coarse-graining increases the approximate state error, here we will attempt to restore positivity without the full-scale coarse-graining of the Redfield equation that we discussed in Sec. II B.
First, we rewrite the Redfield equation (Eq. 11) in terms of the geometric mean spectral density, instead of the arithmetic one. For example, we write (

The Redfield equation becomes
where we introduce the "detuning" function, Note that f (ω, ω) = 0, which is why we call it the detuning function.
For the Ohmic bath at zero temperature, the images in Fig. 1 a) and d) display the geometric mean g(ω, ω ′ ) = γ(ω)γ(ω ′ ) and the detuning function magnitude versus system frequencies, respectively. The entities f and g, that we call kernels, are independent of the system state and the operator A. Without any coarse-graining, the kernels are overall comparable in magnitude, which is unfortunate, because, if the detuning function were negligibly small compared to the geometric mean, we could neglect the second line in Eq. 26, and end up with a Markovian master equation.
Note that the detuning function is large in the frequency quadrants ωω ′ < 0.
In spite of the significant value of the detuning function, we can make a compromise following the steps discussed next. We recast the coarse-grained Eq. 19 into the form of Eq. 26, whereg is the coarse-grained geometric mean, 10, respectively. Gray scale range is the same between the panels, but reduced by factor of two relative to that in panels a)-c).
andf is the coarse-grained detuning function, Let us investigate the effects of coarse-graining on times T 0 = 3/ω c and 10/ω c . The coarsegrained geometric mean spectral density, displayed in Figs. 1 b) and c), is significantly suppressed in the frequency range |ω − ω ′ | > 1/T 0 , but remains mostly unchanged at frequencies near the diagonal ω = ω ′ . The coarse-grained detuning function is similarly suppressed in the frequency range |ω − ω ′ | > 1/T 0 , as shown in Figs. 1 e) and f). In contrast to the geometric mean, the coarsegrained detuning function is strongly suppressed in its entirety because, by design, the detuning function is small at small |ω − ω ′ |. The key here is that the characteristic coarse-graining time that suppresses the detuning function relative to the geometric mean is determined by the bath We use the Frobenius norm to quantify the effect of coarse-graining on the kernels, for some function y(ω, ω ′ ). Fig. 2 a) displays the norms of the coarse-grained kernels. The norm of the coarse-grained detuning function drops much faster than that of the geometric mean. The ratio of the norms of the coarse-grained kernels, as a function of the coarse-graining time is displayed in Figs. 2 b) and c). The time on which the norm ratio drops by 50% is indicated by the arrow, and corresponds to T 0 ≈ 1.75/ω c . At large T 0 , the ratio of the coarse-grained kernels is inversely The best inverse fit is shown by the black line in Fig. 2 c), leading to c N = 0.658. We find similar properties of the norm ratio if we use the trace-norm or the operator norm.
This analysis conveys a picture of dissipative quantum dynamics. The effects of the detuning function average out relative to the effects of the geometric mean density on timescales governed by and somewhat larger than τ c . If we consider system dynamics on timescales much longer than τ c , we can neglect the detuning function.
Next, we introduce the generator matrix, (where √ γ is defined as the matrix with elements √ γ ij ), which leaves us with the Markovian that is more general than the Davies-Lindblad equation. As in the RWA, only positive frequencies ω ij have nonzero elements L ij at zero temperature.
We call this a "Markovian arithmetic geometric mean approximation", in short MAGMA, because complete positivity of master equation is restored, by replacing the arithmetic mean of spectral densities with the geometric mean in the dissipative part of the system dynamics. This shift in the mean changes how dissipation is accounted for at zero temperature. Speaking qualitatively, the arithmetic mean adds dissipative channels in parallel, while the geometric mean is akin to adding them in series. The shift is particularly significant in the frequency quadrants ωω ′ < 0, where at zero temperature, the geometric mean rate is zero, while the arithmetic one can be significant.
In appendix VII we obtain an inequality for the trace distance between the solutions of the Redfield and MAGMA equations, for an arbitrary parameter T 0 . The error bound on the trace distance is obtained by minimizing the RHS versus T 0 , which leads to and is equivalent to the CGSE error bound given by Eq. 25 (after replacing τ c with 1/ω c , and τ sb with τ r as discussed repeatedly in sections II A and II B). We conclude that MAGMA is at least as accurate as the CGSE.
MAGMA can accurately describe the quantum correlations induced in the system by the heat bath. This is primarily due to the unitary component of the system evolution, governed by the renormalized Hamiltonian given by Eq. 9, as we will illustrate in the following sections. The coherences can also by induced by the dissipator L † ρL − {LL † , ρ}/2, but we find that the unitary process dominates. A master equation similar to Eq. 34 was investigated recently. In particular, Eq. 33 is introduced phenomenologically, as an ansatz for quantum kinetics equations. However, the renormalization of the system Hamiltonian was not utilized. MAGMA finds harmony in the renormalization of the system Hamiltonian with the geometric mean approximation. The former changes the loss terms in the master equation, which enables the implementation of the geometric mean to both the gain and loss terms, in a balanced way. This can lead to outstanding agreements with true quantum dynamics, as will be shown in the following sections.

III. 3-level system
In this section, we test the accuracy of the approximation on an exactly solvable model system complex enough to display quantum correlations between the states of the system due to the coupling to the heat bath. For the test, we utilize a V-type 3-level system studied in Ref. 14 .
The unperturbed system eigenenstates are |0 , |2 , and |3 , with the respective eigenenergies E 0 = 0, E 2 , and E 3 . (We skipped the index 1 on purpose.) The total Hamiltonian of the system and the bath of noninteracting linear harmonic oscillators, is where H sb is the system-bath coupling Hamiltonian Here b k are boson annihilation operators, and ω k ≥ 0. The appropriate bath-correlation function for this widths Γ 2 = γ(E 2 ) and Γ 3 = γ(E 3 ). In case A, as shown in Fig. 3 a), the widths are much lower than the level spacing ∆E = E 3 − E 2 , and vice versa in case B [Fig. 3 c)]. where or explicitly where . Lastly, the generator in equation 34 has matrix elements L nm = C nm √ γ nm .
For the initial condition, we suppose that the system is in state |1 . Repeating from Ref. 14 , the time dependent state of the total system in the interaction picture is where |0 b is the bath vacuum. The exact solution follows by numerically solving a system of two Volterra integro-differential equations, where f 2 (t) = exp(iω 20 t)C(t), f 3 (t) = exp(iω 30 t)C(t), f ⋆ g = t 0 f (τ )g(t − τ )dτ . For the system parameters we use g = 0.001, and E 2 = 0.095ω c , and E 3 = 0.105ω c in case A, and E 2 = 0.09975ω c and E 3 = 0.10025ω c in case B. These are chosen to be the same parameters as in Ref. 14 , so that we can directly compare the accuracy of the master equations. Figure 4 displays numerically obtained exact populations of levels 2 and 3 versus time, in parameter cases A and B, respectively. The blue and red curves reproduce previous work. 14 A. Accuracy of the approximate solution Next we investigate the error of the approximate state. We solve Eq. 34 as explained in appendix VIII. The approximate populations of levels 2 and 3 versus time are displayed by the black lines in Fig. 4. The difference between the approximate and the exact solution is barely discernible.
Not shown is the coherence ρ 23 (t), which also agrees well with the exact solution.
. Ref. 14 The trace distance between the approximate and exact solutions is displayed by the black lines in Fig. 5. Also shown by the blue line in Fig. 5(a), is the error of the optimal CGSE, obtained by digitizing the blue line of data in Fig. 3 of Ref. 14 . MAGMA is almost an order of magnitude more accurate than the CGSE, and almost completely captures the population oscillations as seen by the diffuse error with respect to time.
We also determine the trace distance between the solutions of the Redfield and MAGMA equations, and find that the distance is ten times lower than the black lines in Fig. 5. Additionally, solving Eq. 34 without renormalizing the system Hamiltonian (Eq. 42) produces an error much larger than the black lines in Fig. 5. The error increases by a factor of several hundred and demonstrates the necessity to renormalize the Hamiltonian.

B. Results and discussion
Interestingly, Figs. 4 c) and d) show a few initial oscillations that decay rapidly followed by a much slower relaxation process without any oscillations. In case A, as can be seen in Figs. 4 a) and b), the oscillations and populations decay on a similar scale. At first sight, this property is puzzling since the Fermi-golden rule rates are approximately the same in cases A and B. Thinking further, however, the disparity between the relaxation rates is not utterly surprising.
Namely, if levels E 2 and E 3 were degenerate, one could find a basis such that one state has a zero matrix element in the coupling Hamiltonian. For example, if in one basis the coupling Hamiltonian is represented by the matrix in Eq. 39, we can change the basis to |0 , (|2 − |3 )/ √ 2 is the dark state. 14 If the levels are not exactly, but nearly degenerate, the bath will renormalize the states so that one state approaches the dark state resulting in the suppressed relaxation rate.
Here we present a simple procedure that determines the system oscillation frequencies and relaxation rates based entirely on the properties of the renormalized Hamiltonian given by Eq. 43.
The renormalized energies are defined as the nonzero eigenvalues of the matrix in Eq. 43, where E = (E 2 + E 3 )/2, ∆S = S(E 3 ) − S(E 2 ), and the other terms are identical to those in Eq. 43. Let us introduce a parameter dependence of the system Hamiltonian, As a function of λ, the energy levels of H 0 cross at λ = 1, while the renormalized eigenenergies (Eq. 47) have avoided crossing as shown in Fig. 6 a) with the energy gap 2S(E 2 ). The avoided crossing regime is found in the region of energy where the widths γ(E 2,3 ) are comparable to the spacing E 3 − E 2 . In the regime of strong anticrossing, e.g., .
At the center of the crossing, the renormalized states are The last is the dark state.
Here, the avoided crossing is caused by the coupling between the system and the environment of linear harmonic oscillators. We find that not only the states, but also the relaxation rates are strongly affected near such avoided crossings. To obtain the relaxation rates, we consider the initial states of the system ρ(0) = |2 ′ 2 ′ | and |3 ′ 3 ′ | and determine the populations of the renormalized states versus time again by solving Eq. 34. The results are shown in Fig. 7 for case B. In contrast to the time dependence of the unperturbed states shown in Fig. 4, there are no oscillations in Fig. 7.
The goodness of fit to e −Γ ′ t in Fig. 7 indicates that the populations of the renormalized states decay with single effective or renormalized widths.
The width Γ ′ versus λ is displayed in Fig. 6 b). The higher energy state has a strongly suppressed width and becomes dark at λ = 1, while the lower one has enhanced width and increases in brightness by factor of two at λ = 1. Quantum dynamics is fundamentally positive. This suggests that restoring positivity of the Redfield equation may reduce the error relative to the exact solution of the reduced system dynamics The reduction of error, however, is not guaranteed as seen in the RWA. The 3-level quantum system described in the previous section did not have negative Redfield dynamics. In this section, we investigate a 4-level system that is also exactly solvable. The chosen system gives negative Redfield dynamics allowing us to study how restoring positivity affects the error.

A. Setting up the Hamiltonian
We consider a 4 level system described by the Hamiltonian or briefly where B is given by Eq. 40. The Hamiltonian commutes with the operator which greatly simplifies finding the exact quantum dynamics. The levels are labeled in ascending order, e.g., (0,1,2,3). For the energies, we assume E 0 < E 2 < E 3 < E 1 .

B. Exact solution
We seek the solution of the time dependent global wavefunction in the form similar to Eq. 44, Here we assume the initial condition |ψ(0) = |2 ⊗ |0 b . Following the steps presented in Ref. 14 , we arrive to the pair of Volterra integro-differential equations (Eqs. 45 and 46). Since the derivation of these equations mostly repeats that of the 3-level system, we omit it here. The 4-level system differs from the 3-level system in the following expressions: f 2 (t) = [exp(iω 20 t) +

exp(iω 21 t)]C(t) and f 3 (t) = [exp(iω 30 t) + exp(iω 31 t)]C(t).
The solutions d k (t) and e k (t) can be obtained from the Schrödinger equation in the interaction picture, i|ψ(t) = H sb (t)|ψ(t) , which can formally be integrated to |ψ(t) = |ψ(0) − i t 0 dτ H sb (τ )|ψ(τ ) , and leads to and The reduced density matrix of the system, ̺ = T r b |ψ(t) ψ(t)|, has block diagonal form where r 00 = k |d k | 2 , r 01 = k d k e ⋆ k , and r 11 = k |e k | 2 . These sums can be evaluated in terms of the bath-correlation function. For example, where the sum on the RHS is the bath-correlation function.
The r-block on the upper left side of the reduced state follows from the c-block as We solve these integrals numerically and confirm the accuracy by checking that T r̺ = 1.

C. Master equations
Regardless of what master equation we use, we find that, for the 4-level system, the density matrix decomposes into a direct sum as in Eq. 60. The derivative of the r-block is only a function of the c-block, while the derivative of the c-block does not depend on the r-block. Quantum dynamics within the c-block is very similar to the dynamics of the 3-level system, with excellent agreements between the exact and two approximate solutions. As earlier, the c-block has positive Redfield dynamics.
In contrast to the c-block, different approximations produce very different dynamics in the r-block. Fig. 8 a)  The two-by-two r-block will be positive, if and only if |r 01 | ≤ √ r 00 r 11 .
This means that for a given coherence r 01 , of a positive state, the population r 11 cannot be arbitrarily low. The Redfield equation fails in this respect, since it severely under-populates the highest energy state 1 (Fig. 8 a). The negative eigenvalues of the state it produces, are displayed in Fig. 8 c). Black circles in Fig. 8 b) correspond to the lowest three minima in Fig. 8 c). We notice that the state is the most negative, approximately when the coherences produced by the Redfield equation have maximum magnitude, and vice versa.
In contrast to the Redfield equation, MAGMA has trivially balanced populations and coherences. They are both zero, and therefore satisfy the inequality of Eq. 65. The balance is restored by the shift from the arithmetic to geometric mean spectral density.
Next we discuss the accuracy of the dynamics described by the two master equations. Fig. 8 d) displays the trace distances between the r-blocks of the exact and approximate solutions showing that MAGMA has smaller error than the Redfield equation. The reason is that the latter has an erroneous coherence. The distance between the exact coherence and zero is smaller, by roughly a factor of √ 2 than that between the exact and the erroneous coherence of similar magnitude.
In summary, for this very simple 4-level system, MAGMA performs better than the Redfield master equation by maintaining complete positivity and producing a smaller error. In the next section, we examine a much more complex many body system with very small level spacing, that cannot be solved exactly, and study the distances between the Redfield and MAGMA states.

V. Magnetic relaxation in a ferromagnetic spin-chain
Here we consider a spin-chain of 25 spin-1/2 particles to study the quantum dynamics in large open many-body systems. Each particle couples to its nearest neighbor via the ferromagnetic exchange interaction and to each other through the long-range, magnetic dipole-dipole interaction. This system has a rapidly increasing density of states with energy and is rich in many-body phenomena, including, the energy gap (ferromagnetic resonance), collective excitations (standing spin-waves), and topological excitations (domain walls). There is a corresponding plethora of relaxation processes, including Gilbert damping, emission of spin-waves, and domain wall motion, spanning a wide range of relaxation times, that can be exponential in system size. The focus here will be on the relaxation of the density matrix in response to a macroscopic displacement of the magnetization, which is initially set perpendicular to the chain (easy) axis. We investigate the difference between the Redfield, RWA, and MAGMA master equations, at the early stage of the relaxation cascade, before the slow moving domain walls nucleate.

A. Magnetic Hamiltonian
The system is modeled by the Heisenberg model with dipole-dipole coupling, The Hamiltonian commutes with S z = i S iz , but does not commute with S 2 , where S = i S i .
Here n = 25 and the eigenvalues display Kramers degeneracy. In addition, the Hamiltonian is For parameters we choose J = 400 and ǫ d = 6. With no anisotropy, (e.g., ǫ d = 0), the spacing between the lowest lying spin-multiplets is δ = 2.5. Each multiplet represents a different standing spin wave, while the population of the standing wave varies within the corresponding multiplet.
In the case ǫ d = 6, the ground state is that of a uniaxial ferromagnet, |S, S z = ±S , (chain axis is along z-direction), the ground state energy is E g = −2485.28, and there is an anisotropy gap ∆ = 20 in the excitation spectrum, as can be seen in Fig. 9. The gap is the ferromagnetic resonance. We work in the regime where ∆ > δ, which is the regime common in devices. We are interested in the low energy collective magnetic dynamics. It suffices to find the eigenvalues and eigenstates, spanning the energy range between the ground state and the Néel barrier, as we will see shortly. In particular, let us consider the initial state with the magnetization fully oriented along the x-direction. To obtain this initial state, we truncate S x to the vector space on the lowest 1728 eigenstates and find the eigenstate of S x with the maximum eigenvalue. In that state, we obtain S x ≈ 12.41, which is lower than 12.5 because of the truncation error. Since this value is very close to 12.5, it shows that the vector space is sufficiently large to describe uniformly magnetized states, and the states accessible by relaxation from those states. The expectation value of the energy in the initial state is approximately the same as the energy at the top of the blue parabola in Fig. 9, while the energy uncertainty is rms(H 0˙) = 11.6.
The energy levels under the parabola in Fig. 9 are noteworthy. They are grouped into branches, which have weak dependence of energy with S z . On these branches, energies have two-fold near degeneracy (not visible in Fig. 9) between even and odd spin states. At small |S z |, these states are very similar to domain walls located at distance S z from the middle of the chain. The energy of the domain wall doesn't change much if the domain moves near the middle of the chain, explaining the flatness of the branches. At this time, it would be a distraction to delve into these states. What matters at this point is that there is an energy gap above the ground state, a major increase in the density of states at the top of the barrier, and a wide distribution of relaxation times due to large variety of the excitations making this many-body system a good example to investigate complex open quantum dynamics.

B. Coupling to the heat bath
Each spin is coupled to a noisy magnetic field, analogous to the environment in the celebrated Brown model. 24 In our case, there are 75 independent baths of linear harmonic oscillators, 3 for each spin, one for each direction of the field. The coupling Hamiltonian is Here are Hermitian bath operators, and we assume slow bath at zero temperature with Ohmic spectral density and cut-off frequency ω c = 6∆.

The Redfield equation for the spin-chain is
where S iα,f = S iα • F , and F is given after Eq. 5.

The corresponding equation in MAGMA is
with the renormalized Hamiltonian, The Lindblad generators are obtained from Eq. 33, L iα = S iα • √ γ. We solve these equations by numeric iterations, in the truncated Hilbert space, as explained in appendix VIII. For completeness, we also find the solutions using the RWA.

C. Relaxation Processes
First we study the relaxation of the unform spin state, with the initial spin oriented perpendicular to the chain. We assume g = 0.000133 for the spectral density in each bath. Since there are 75 independent baths, the total coupling is g tot = 75g = 0.01. In Fig. 9 b) the bottom axis is time in units of T f m , where T f m = 2π/∆ is the period of the ferromagnetic resonance. The figure displays rapid decay of the perpendicular magnetization, which in this case is mostly due to the unitary evolution, caused by the energy uncertainty of the initial state [e.g., T rρ 2 is not much ) " lower than one in Fig. 9 b)]. The energy of the system H 0 is shown in Fig. 9 c). It decays due to the dissipative coupling to the heat bath. The small drop of H 0 , relative to the initial energy, indicates that the figure displays an early stage of the relaxation cascade.
Next we find the time scale τ r , on which the heat bath affects the state. We rotate into the interaction picture ̺ = exp(iH 0 t)ρ exp(−iH 0 t), and determine τ r as this, where the subscript R indicates the state is obtained by solving the Redfield equation. From now, we approximate this as since ||̺ R || 1 is only slightly larger than one, due to the negativity of the Redfield state, in the examples we study in this paper. Fig. 10 a) displays 1/τ r versus time, for g tot = 0.005 and 0.01.
The state relaxation rate 1/τ r scales with the g tot , as expected. We find the rate 1/τ r by timeaveraging, which calibrates the rate according to 1/τ r = 5.91g tot .
The trace distance between the solutions of the Redfield and MAGMA is displayed in Fig. 10 b). MAGMA produces states that are much closer to the Redfield states, compared to that of the RWA. To get a better sense of this distance, in Fig. 11 a) and b) we display the trace distance versus time and the coupling strength g tot . On panel b) we scale the time with τ r . Next, in Fig. 11 c), we show the maximum trace distance between the solutions of the Redfield and MAGMA equations versus the state relaxation rate 1/τ r . The figure displays asymptotic linear scaling between the trace distance and the state relaxation rate at low 1/τ r . The best linear fit on the lowest four points is 1 2 ||ρ Redf ield − ρ magma || 1,max ≈ 6.41 1 ω c τ r , and is shown by the blue line. At large g tot , as the trace distance begins to approach the saturation value of one, the dependence becomes sublinear.
The linear scaling between the maximum trace distance and the state rate τ −1 r suggests that the error bound given by Eq. 36 is not tight. That is, it seems that the high accuracy of MAGMA is not due to a smaller prefactor in the expression for the error of the coarse-grained stochastic master equation. Rather, the linearity of the trace distance with the state relaxation rate with the slope not far off the Born-Markov error, which is also linear with the state relaxation rate, suggests that MAGMA has the error set by the Born-Markov error (Eq. 7).

VI. Conclusion
In this article we derive a Markovian quantum master equation from the Redfield master equation by truncating some, but not all, terms that average out over a time-scale typical of the bath.
We estimate that the state error resulting from this approximation is smaller than or comparable to the error in the coarse-grained stochastic master equation. However, in practice we find the error to be much lower than the latter and close to the Born-Markov error. The Markovian master equation is as simple to use as the Redfield equation. In one example where the exact quantum dynamics of the open quantum system is known, we find that the error of the Markovian master equation is lower than the error of the Redfield equation.
The key to achieving high accuracy is in the renormalization of the system Hamiltonian due to the coupling to the heat bath. This renormalization must be done before the arithmetic to geometric mean spectral density approximation, or MAGMA, to ensure harmony between the unitary and dissipative quantum dynamics.
There is no fundamental reason that a completely positive master equation cannot have an error as small as the Born-Markov error. The examples that we investigate in this paper clearly demonstrate that. The next step will be to extend the equation to time-dependent Hamiltonian, which is important to conceive future experiments for quantum information and nonequilibrium quantum dynamics in many-body systems. We thank Jason Dark, Elyana Crowder, Brian Kennedy, and Glen Evenbly for discussions. This research was supported by DOE contract DE-FG02-06ER46281, including development of the numerical simulations of quantum dynamics in magnets. 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 systems.

VII. Appendix: Error Bound
Here we estimate the error bound for the trace distance between the solutions of the Redfield and MAGMA equations, used in Sec. II D. This estimate falls short of a full proof but is meant to be a plausibility argument. Denote by L{ρ} and F {ρ} the superoperators on the first and second lines on the RHS of Eq. 10. L is the Lindbladian, and we refer to F as the detuner. The Redfield and MAGMA equations are The respective coarse-grained equations are Here L =ĜL, F =ĜF , andĜ is defined in Eq. 15. According to Lemma 1 in Ref. 16 , the error between the solutions of Eqs. 76 and 77, at a fixed time on the order of the state relaxation time, The trace distance between the solutions of the Redfield and MAGMA master equations can be estimated using the triangle inequality, Using Eqs. 22 and 78, this leads to the trace distance Here we assumed that τ r in the Redfield and MAGMA equations is approximately the same, as confirmed in all examples studied in this paper. It is also sensible to assume that, because otherwise MAGMA would not be a good approximation. We make the last assumption that the norms of the superoperators and the kernels linearly scale, where the terms on the RHS are defined in Eqs. 31, 30 and 29. Then, using the asymptotic dependence of the norm ratio given by Eq. 32, we arrive at Eq. 35.

VIII. Appendix: Numerical Iteration of Master Equations
We discuss how we solve a master equation within an arbitrary error (ǫ) in norm distance, given the initial condition ρ(0). Here L is a time independent superoperator. The algorithm simulates the power series expansion of e Lt , e.g., as follows. Given the state ρ(t) at time t, and time increment dt, let ρ d = ρ 0 = ρ(t). Start the loop m = 1, 2, ... ,. For each m, let ρ d → L{ρ d }, followed by ρ 0 → ρ 0 + (dt m )ρ d /m!. Finish the loop when ||(dt m )ρ d /m! || < ǫ. All results discussed in this article have accuracy of at least three significant digits.