Accuracy Assessment of Perturbative Master Equations -- Embracing Non-Positivity

TL;DR

This paper evaluates the accuracy of various perturbative master equations for open quantum systems, highlighting that non-positivity violations can serve as indicators of the breakdown of weak coupling assumptions.

Contribution

It provides a comprehensive assessment of multiple master equations, emphasizing the significance of non-positivity as a validity indicator in quantum dynamics modeling.

Findings

Redfield Equation with time-dependent coefficients is most accurate.

Positivity violations occur only when perturbative methods become invalid.

Non-positivity should be viewed as a sign of the breakdown of weak coupling assumptions.

Abstract

The reduced dynamics of an open quantum system obtained from an underlying microscopic Hamiltonian can in general only approximately be described by a time local master equation. The quality of that approximation depends primarily on the coupling strength and the structure of the environment. Various such master equations have been proposed with different aims. Choosing the most suitable one for a specific system is not straight forward. By focusing on the accuracy of the reduced dynamics we provide a thorough assessment for a selection of methods (Redfield Equation, Quantum Optical Master Equation, Coarse-Grained Master Equation, a related dynamical map approach and a partial-secular approximation). Whether or not an approach guarantees positivity we consider secondary, here. We use two qubits coupled to a Lorentzian environment in a spin-boson like fashion modeling a generic situation with various system and bath time scales. We see that, independent of the initial state, the simple Redfield Equation with time dependent coefficients is significantly more accurate than all other methods under consideration. We emphasize that positivity violation in the Redfield formalism becomes relevant only in a regime where any of the perturbative master equations considered here are rendered invalid anyway. This implies that the loss of positivity should in fact be welcomed as an important feature: it indicates the breakdown of the weak coupling assumption. In addition we present the various approaches in a self-contained way and use the behavior of their errors to provide further insight into the range of validity of each method.

Equations

$$\begin{array}{c}H=H_{\mathrm{sys}}+L\otimes \sum _{\lambda }{g}_{\lambda }(a_{\lambda }+a_{\lambda }^{†})+\sum _{\lambda }{\omega }_{\lambda }{a}_{\lambda }^{†}{a}_{\lambda }\\ H_{\mathrm{sys}}=\frac{{\omega }_A}{2}{\sigma }_x^A+\frac{{\omega }_B}{2}{\sigma }_x^{B}L=\frac{1}{2}({\sigma }_z^A+{\sigma }_z^B).\end{array}$$

$$J(w)=\eta \frac{\gamma }{{\gamma }^2+({\omega }_0-\omega {)}^2}.$$

$$\alpha (\tau )=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega J(\omega )e^{-\mathrm{i}\omega \tau }=\eta e^{-\gamma |\tau |-\mathrm{i}{\omega }_0\tau },$$

$$F(\omega )$$

$$S(\omega )$$

$$\begin{gathered} H=H_{\mathrm{sys}}+\sqrt{\eta }L(a^{†}+a)+{\omega }_0{a}^{†}a \\ +\sum _{\lambda }{c}_{\lambda }(a^{†}{b}_{\lambda }+a{b}_{\lambda }^{†})+\sum _{\lambda }{\omega }_{\lambda }{b}_{\lambda }^{†}{b}_{\lambda } \end{gathered}$$

$$\dot{P}$$

$$H_{\mathrm{sys}}^P$$

$$\dot{\tilde{\rho }}(t)=-{\int }_0^t\mathrm{d}s{\mathrm{Tr}}_{\mathrm{env}}[\tilde{V}(t),[\tilde{V}(s),\tilde{\rho }(s)\otimes {\rho }_{\mathrm{env}}]].$$

$$\dot{\tilde{\rho }}(t)=-{\int }_0^t\mathrm{d}s\left(\alpha (t-s)[\tilde{L}(t),\tilde{L}(s)\tilde{\rho }(s)]+h.c.\right)$$

$$\alpha (\tau )$$

$$\displaystyle=\sum_{\lambda}g_{\lambda}^{2}\big{(}(2\bar{n}(\beta\omega_{\lambda})+1)\cos(\omega_{\lambda}\tau)-\mathrm{i}\sin(\omega_{\lambda}\tau)\big{)}\quad.$$

$$\begin{gathered} \dot{\rho }(t)=-\mathrm{i}[H_{\mathrm{sys}},\rho (t)] \\ +{\int }_0^t\mathrm{d}\tau \left(\alpha (\tau )[\tilde{L}(-\tau )\rho (t),L]+h.c.\right). \end{gathered}$$

$$L_{\omega }=\sum _{ϵ,{ϵ}^{\prime }:{ϵ}^{\prime }-ϵ=\omega }|ϵ⟩⟨ϵ|L|{ϵ}^{\prime }⟩⟨{ϵ}^{\prime }|$$

$$L(t)=e^{\mathrm{i}{H}_{\mathrm{sys}}t}L{e}^{-\mathrm{i}{H}_{\mathrm{sys}}t}=\sum _{\omega ={ϵ}^{\prime }-ϵ}{e}^{-\mathrm{i}\omega t}{L}_{\omega }.$$

$$\begin{gathered}\dot{\rho}(t)=-\mathrm{i}[H_{\mathrm{sys}},\rho(t)]+\sum_{\omega}\Big{(}F(\omega,t)[L_{\omega}\rho(t),L]+\mathrm{h.c.}\Big{)}\\ F(\omega,t)=\int_{0}^{t}\mathrm{d}\tau\,\alpha(\tau)e^{\mathrm{i}\omega\tau}\quad.\end{gathered}$$

$$F(\omega ,t)=\eta \frac{\gamma +\mathrm{i}(\omega -{\omega }_0)}{{\gamma }^2+({\omega }_0-\omega {)}^2}\left(1-e^{-(\gamma +\mathrm{i}({\omega }_0-\omega ))t}\right).$$

$$\dot{\tilde{\rho }}(t)=\sum _{\omega ,{\omega }^{\prime }}{e}^{-\mathrm{i}(\omega -{\omega }^{\prime })t}F({\omega }^{\prime },t)[L_{{\omega }^{\prime }}\tilde{\rho }(t),L_{\omega }^{†}]+h.c..$$

$$\dot{\rho}(t)=-\mathrm{i}[H_{\mathrm{sys}}+\sum_{\omega}S(\omega)L^{\dagger}_{\omega}L_{\omega},\rho(t)]\\ +\sum_{\omega}\Big{(}J(\omega)[L_{\omega}\rho(t),L_{\omega}^{\dagger}]+\mathrm{h.c.}\Big{)}\quad.$$

$$L_{{\omega }_{A/B}}=\frac{1}{2}|{\psi }_{-}{⟩}^{A/B}⟨{\psi }_{+}{|}^{A/B}=L_{-{\omega }_{A/B}}^{†}$$

$$L_{\omega }=\frac{1}{2}\left(|{\psi }_{-}^A⟩⟨{\psi }_{+}^A|+|{\psi }_{-}^B⟩⟨{\psi }_{+}^B|\right)=L_{-\omega }^{†}.$$

$$\dot{\tilde{\rho }}(t)=\sum _{\bar{\omega },{\bar{\omega }}^{\prime }}{e}^{-\mathrm{i}(\bar{\omega }-{\bar{\omega }}^{\prime })t}F({\bar{\omega }}^{\prime },t)[L_{{\bar{\omega }}^{\prime }}\tilde{\rho }(t),L_{\bar{\omega }}^{†}]+h.c.$$

$$\dot{\rho}(t)=-\mathrm{i}[H_{\mathrm{sys}}+\sum_{\bar{\omega}}S({\bar{\omega}})L^{\dagger}_{\bar{\omega}}L_{\bar{\omega}},\rho(t)]\\ +\sum_{\bar{\omega}}\Big{(}J({\bar{\omega}})[L_{\bar{\omega}}\rho(t),L_{\bar{\omega}}^{\dagger}]+\mathrm{h.c.}\Big{)}\quad.$$

$$\begin{gathered} {\tilde{\rho }}_{\mathrm{tot}}(t+\tau )-{\tilde{\rho }}_{\mathrm{tot}}(t)\approx -\mathrm{i}{\int }_t^{t+\tau }\mathrm{d}s[\tilde{H}(s),{\tilde{\rho }}_{\mathrm{tot}}(t)] \\ -{\int }_t^{t+\tau }\mathrm{d}s{\int }_t^s\mathrm{d}u[\tilde{H}(s),[\tilde{H}(u),{\tilde{\rho }}_{\mathrm{tot}}(t)]], \end{gathered}$$

$$\tilde{\rho}(t+\tau)-\tilde{\rho}(t)\approx-\int_{t}^{t+\tau}\mathrm{d}s\,\int_{t}^{s}\mathrm{d}u\,\Big{(}\alpha(s-u)\\ \times[\tilde{L}(s),\tilde{L}(u)\tilde{\rho}(t)]+\mathrm{h.c.}\Big{)}=:\mathcal{Z}_{\tau}\tilde{\rho}(t)$$

$$\dot{\tilde{\rho }}(t)\approx \frac{\tilde{\rho }(t+\tau )-\tilde{\rho }(t)}{\tau }\approx \frac{{\mathcal{Z}}_{\tau }}{\tau }\tilde{\rho }(t).$$

$$\dot{\tilde{\rho}}(t)=\frac{\mathcal{Z}_{\tau}}{\tau}\tilde{\rho}(t)=-\frac{1}{\tau}\sum_{\omega,\omega^{\prime}}\Big{(}e^{\mathrm{i}(\omega-\omega^{\prime})t}G(\omega,\omega^{\prime},\tau)\\ \times[L_{\omega^{\prime}},L_{\omega}^{\dagger}\tilde{\rho}(t)]+\mathrm{h.c.}\Big{)}$$

$$G(\omega ,{\omega }^{\prime },\tau )={\int }_0^{\tau }\mathrm{d}s{\int }_0^s\mathrm{d}u\alpha (s-u)e^{\mathrm{i}({\omega }^{\prime }s-\omega u)}$$

$$\begin{gathered} G(\omega ,\omega ,\tau )=\frac{\eta }{\gamma +\mathrm{i}({\omega }_0-\omega )}\tau \\ +\frac{\eta }{(\gamma +\mathrm{i}({\omega }_0-\omega ){)}^2}\left(e^{-(\gamma +\mathrm{i}({\omega }_0-\omega )\tau }-1\right) \end{gathered}$$

$$G(\omega,\omega^{\prime},\tau)=\frac{\eta}{\gamma+\mathrm{i}(\omega_{0}-\omega)}\Bigg{[}\frac{\mathrm{i}}{\omega-\omega^{\prime}}\left(e^{-\mathrm{i}(\omega-\omega^{\prime})\tau}-1\right)\\ +\frac{1}{(\gamma+\mathrm{i}(\omega_{0}-\omega^{\prime}))}\left(e^{-(\gamma+\mathrm{i}(\omega_{0}-\omega^{\prime})\tau}-1\right)\Bigg{]}\quad.$$

$$\underset{\tau \to \infty }{\lim }\frac{G(\omega ,{\omega }^{\prime },\tau )}{\tau }=(J(\omega )+\mathrm{i}S(\omega )){\delta }_{\omega ,{\omega }^{\prime }}.$$

$$\tilde{\rho }(t)=(\text{mathds}1+{\mathcal{Z}}_t)\tilde{\rho }(0),$$

Full text

I Introduction

The non-unitary dynamics of an open quantum system is of great interest for many fields in physics and chemistry, where dissipation and decoherence have to be considered. Solving the microscopic model of the whole, system plus environment, with regard to the exact reduced dynamics is in general still a difficult task (see, e.g. Makri (1995); Thorwart and Jung (1997); Beck et al. (2000); Wang and Thoss (2003); Ishizaki and Tanimura (2005); Tanimura (2006); Suess et al. (2014); Hartmann and Strunz (2017)). However, in the weak coupling regime time local master equations can be derived from the microscopic model resulting in approximate solutions for the reduced dynamics Cohen-Tannoudji et al. (1998); Breuer and Petruccione (2007); Weiss (2008). Such master equations can easily be solved with standard numerical methods. However, that advantage is to some extent dissolved by the lack of a possibility to estimate the error of the approximations from within the method. Consequently, many variants of microscopically motivated master equations have been proposed pending a rigorous analysis of their applicability. In the work presented here we provide that missing study. By considering a Lorentzian environment we are able to obtain the exact reduced dynamics by means of the pseudo-mode method. Henceforth, for the exemplary system of two spins (qubits) coupled in spin-boson like fashion Leggett et al. (1987) to the environment we can explicitly calculate the error of the reduced dynamics obtained from the master equations. To assess each master equation by its accuracy, we calculate an initial state independent error bound and plot its behavior as a function of bath correlation time and coupling strength, the two quantities primarily influencing the approximations involved.

In the following we shortly review the motivation for a selection of master equations. The least approximations are required to derive the Redfield Equation Redfield (1957); Breuer and Petruccione (2007); Weiss (2008). Since that equation is not of Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) form, positivity of the reduced state is in general not guaranteed Suárez et al. (1992); Pechukas (1994); unphysical negative eigenvalues of the density operator may occur after some time. However, this shortcoming does not imply that the dynamics obtained from the Redfield Equation is of little use. In particular, sufficiently weak coupling and a fast decaying bath correlation function (BCF) justify the approximations made, which do render the solution of the Redfield Equation valid within a certain error range. From a practical point of view, this motivates the wide application of the Redfield Equation and its variants (see for example Refs. Kohen et al. (1997); Kondov et al. (2001); Egorova et al. (2003); Nitzan (2006); Schröder et al. (2007); Timm (2008); Montoya-Castillo et al. (2015); Bricker et al. (2018) for recent quantum chemical, condensed matter and quantum optics applications).

From a more conceptual open quantum system point of view, the lack of (complete) positivity implied by the Redfield Equation results in a rejection of the method Benatti and Floreanini (2005); Rivas et al. (2010). However, positivity preservation can be enforced by further approximations.

The most prominent and seasoned additional approximation is the so-called rotating-wave or secular approximation (RWA) Lindblad (1976); Lidar et al. (2001); Breuer and Petruccione (2007); Weiss (2008). The resulting Quantum Optical Master Equation is of GKSL-type and therefore ensures completely positive dynamics. The applicability of the RWA, however, requires a sufficiently weak coupling of the system to the environment such that the dynamics of the reduced state in the interacting picture takes place on a much larger timescale than the timescale set by the differences of the transition frequencies of the system Hamiltonian.

If, for example, the RWA is justified for a single qubit ($H_{\mathrm{sys}}=\frac{1}{2}\omega_{A}\sigma_{x}^{A}$) coupled to an environment, it might not be the case for two slightly detuned qubits ($H_{\mathrm{sys}}=\frac{1}{2}\omega_{A}\sigma_{x}^{A}\mathds{1}^{B}+\frac{1}{2}\omega_{B}\mathds{1}^{A}\sigma_{x}^{B}$) coupled simultaneously to the same environment111Throughout this paper we use units where $\hbar=k_{\mathrm{B}}=1$.. The Hamiltonian of this four dimensional system yields a frequency difference of the order of the detuning $\omega_{A}-\omega_{B}$ which introduces a new, presumably larger timescale, requiring an even weaker coupling to the environment for the RWA to be applicable. Recent studies Benatti et al. (2010); Ma et al. (2012); Eastham et al. (2016); Dodin et al. (2018) address this topic and reassure that the RWA yields significant deviation from the exact dynamics. When increasing the dimension of the system Hilbert space this problem is prone to become even more significant.

Consequently, a master equation of GKSL-form which can be derived from the microscopic model without RWA seems desirable. A coarse-graining approach leading to the so-called Coarse-Grained Master Equation Schaller and Brandes (2008); Benatti et al. (2009, 2010); Majenz et al. (2013), with a coarse-graining time $\tau$ as a free parameter, may accomplish that task. It requires no direct RWA which suggests that this method could outperform the usual Quantum Optical Master Equation. Moreover, a completely positive map, here called ExpZ Map, can be constructed from the $\tau$-dependent generator of the Coarse-Grained Master Equation, which yields the correct dynamics for short times, while recovering the long time dynamics of the usual Quantum Optical Master Equation Majenz et al. (2013); Rivas (2017). In addition, the RWA may be applied only partially Vogt et al. (2013); Jeske et al. (2015); Tscherbul and Brumer (2015) which still yields a master equations of GKSL kind. This approach requires that the relevant transition frequencies cluster such that the spectral density (SD) can be assumed constant for each cluster and that the RWA may be applicable on the level of the cluster frequencies.

Naturally, the question for the best method arises. Although the literature on master equations is vast, an answer based on the deviation from the exact reduced dynamics is missing. As a central part of the current work we address that question by ranking the various methods based on rigorous and initial state independent error bounds for the dynamics of the reduced state. Whether or not the master equation is of GKSL kind is secondary. Scanning a wide range of the environmental parameters coupling strength and bath correlation time, and considering the two-spin-boson model which involves more than a single transition frequency for the system, allows us to draw conclusions about the general applicability of each method.

The extension of the spin-boson model to two qubits is not only interesting from a theoretical point of view by challenging the applicability of the RWA, but also relevant for quantum technologies as it serves as a basic building bloc to implement quantum information tasks Imamoğlu et al. (1999); Clarke and Wilhelm (2008); Ladd et al. (2010).

Ranking the various methods based on the error bound shows that the Redfield Equation performs best. This is in line with the fact that the Redfield Equation involves the least approximations. When using the Redfield Equation with time dependent coefficients, we affirm that positivity problems of the reduced dynamics do only occur in a regime where the weak coupling formalism is not applicable anyway. There is, thus, no reason to abandon the Redfield Equation on the basis of positivity considerations. We find that ranking the other methods considered here, all being of GKSL type, in a general way is impossible since their quality depends in a more sophisticated way on the environmental parameters and/or the properties of the system.

Further, we examine the different scaling of the error with the coupling strength. In addition, we elucidate how the Quantum Optical Master Equation differs in the resonant case from the more general detuned case and show that in the detuned case the correlations between the two qubits are strongly effected by the RWA. We find that the ExpZ Map barely resolves the deficiency of the RWA, whereas the Coarse-Grained Master Equation does so in the relatively strong coupling and short bath correlation time regime. The partial RWA, which is of the structure of the Quantum Optical Master Equation for resonant qubits, seems the best candidate to deal with the detuned case while guaranteeing positivity. However, one should bear in mind that the partial RWA approach relies on the specific energy spectrum of the system Hamiltonian whereas the Coarse-Grained Master Equation is independent of that.

The manuscript is structured as follows. In Sec. II we briefly introduce the notation for the two-spin-boson model and present its solution in Sec. III. The exact solution in terms of a single pseudo-mode is explained, followed by the various approximative master equations. The results in Sec. IV begin with a general ranking of the methods based on an initial state independent error bound. A discussion of the positivity of the reduced dynamics for the Redfield Equation follows, where the advantage of the time dependent coefficients over the Redfield Equation with asymptotic rates is highlighted. Next, the linear scaling of the error is shown briefly and the influence of the coarse-graining parameter $\tau$ is discussed. Finally, the particular effect of the various master equations on the delicate correlation dynamics within the two qubits is shown. We close with a summary and conclusions.

II Two-Spin-Boson Model

The Hamiltonian for two qubits coupled to the same environment, which we will refer to as two-spin-boson model, takes the usual form for an open quantum system with a collective Hermitian operator $L=L^{\dagger}$, coupling the two spins to a common bath of harmonic oscillators,

[TABLE]

The coupling constants $g_{\lambda}$ and the oscillator frequencies $\omega_{\lambda}$ define the SD $J(\omega)=\pi\sum_{\lambda}g_{\lambda}^{2}\delta(\omega-\omega_{\lambda})$. For the continuous environment we choose a single Lorentzian-like SD 222At first glance it seems unphysical to include negative frequencies. However, the physical meaning can be restored when viewing the Lorentzian SD as mathematical vehicle to conveniently model a non-zero temperature BCF with microscopically defined SD $J_{0}$: $\alpha_{\beta}(\tau):=\frac{1}{\pi}\int_{0}^{\infty}\mathrm{d}\omega\,J_{0}(\omega)[\coth(\beta\omega/2)\cos(\omega\tau)-\mathrm{i}\sin(\omega\tau)]\rightarrow\frac{1}{\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega\,J(\omega)e^{-\mathrm{i}\omega\tau}$

with central frequency $\omega_{0}$, width $\gamma$ and overall coupling strength $\eta$

[TABLE]

The corresponding BCF takes the very pleasant form of an exponential

[TABLE]

which allows to easily calculate its half-sided Fourier transform

[TABLE]

a function occurring in various master equation approaches.

III Exact Numerical Solution and Master Equations

III.1 Exact Solution - Pseudo Mode

In terms of the reduced dynamics, the general open quantum system Hamiltonian (Eq. (1)) with Lorentzian SD is equivalent to a pseudo-mode model, where the system couples to a single harmonic oscillator with frequency $\omega_{0}$ which in turn is coupled with coupling strength $\sim\gamma$ to an environment with a flat SD Imamoglu (1994); Garraway (1997). In other words, the Hamiltonian

[TABLE]

with $J_{b}(\omega)=\pi\sum_{\lambda}|c_{\lambda}|^{2}\delta(\omega-\omega_{\lambda})=\gamma$ leads to the same reduced dynamics for the system part.

As a consequence of the flat SD for the $b$-modes the imaginary part of the corresponding BCF vanishes and the real part becomes delta-like: $\alpha_{b}(\tau)=2\gamma\delta(\tau)$. In that case the following master equation of GKSL-type for the state $P$ of the system plus pseudo mode $a$ is known to be exact (see also Ref. Mazzola et al. (2009)).

[TABLE]

Solving the above equation numerically and tracing out the $a$-mode yields the exact reduced state of the two-qubit system $\rho_{\mathrm{sys}}(t)=\mathrm{Tr}_{a}P(t)$. For the following it serves as reference, when comparing the accuracy of the various perturbative master equations.

III.2 Master Equations

The goal of the following master equations is to provide an evolution equation for the reduced state of the open quantum system as given in Eq. (1) for an arbitrary SD. Besides sketching derivations, we also examine the effect of the RWA, distinguishing the resonant and detuned two-qubit case. Further, we discuss the implications of the approximations used by the coarse-graining scheme.

III.2.1 Redfield Equation

To derive the evolution equation for the reduced state Breuer and Petruccione (2007); Kryszewski and Czechowska-Kryszk (2008); Whitney (2008) the Nakajima-Zwanzig projection formalism Nakajima (1958); Zwanzig (1960); Grabert (2006) may be used as starting point. In lowest order of the coupling strength the following expression is obtained

[TABLE]

Here $\tilde{\rho}$ and $\tilde{V}$ denote the reduced state and the interaction Hamiltonian in the interaction picture. Also an initial product state of the form $\rho_{\mathrm{tot}}(0)=\rho(0)\otimes\rho_{\mathrm{env}}$ has been assumed. For the microscopic Hamiltonian Eq. (1) we find explicitly $\tilde{V}(t)=\tilde{L}(t)\otimes F(t)$ with the force operator $F(t)=\sum_{\lambda}g_{\lambda}\left(a_{\lambda}e^{-\mathrm{i}\omega_{\lambda}t}+a^{\dagger}_{\lambda}e^{\mathrm{i}\omega_{\lambda}t}\right)$. Assuming a thermal initial state $\rho_{\mathrm{env}}\sim e^{-\beta H_{\mathrm{env}}}$, the evolution equation Eq. (7) becomes

[TABLE]

with the BCF

[TABLE]

For a BCF decaying faster than the dynamical time scale of the reduced state in interaction picture, $\tilde{\rho}(s)$ may well be approximated by $\tilde{\rho}(t)$ under the integral. Finally, substituting $\tau=t-s$ and transforming back to the Schrödinger picture yields

[TABLE]

The remaining interaction picture contribution can be made explicit by trivially rewriting the coupling operator $L$ in terms of eigenbase projectors of the system Hamiltonian $L=\sum_{\epsilon}|{\epsilon}\rangle\!\langle\epsilon|L\sum_{\epsilon^{\prime}}|{\epsilon^{\prime}}\rangle\!\langle\epsilon^{\prime}|$Breuer and Petruccione (2007). Grouping all terms for a fixed $\omega=\epsilon^{\prime}-\epsilon$ defines

[TABLE]

and allows for the decomposition $L=\sum_{\omega}L_{\omega}$ where $\omega$ runs over all possible transition frequencies of $H_{\mathrm{sys}}$. Consequently, for an operator $L$ in the interacting picture we can write

[TABLE]

Finally, we arrive at the Redfield Master Equation with time-dependent coefficients (Redfield Equation (tdc)):

[TABLE]

For the model BCF as given in Eq. (3), a single exponential, the time dependent coefficients can be evaluated explicitly,

[TABLE]

For a sufficiently fast decaying BCF the asymptotic values $F(\omega,t\rightarrow\infty)=J(\omega)+\mathrm{i}S(\omega)$ may be used instead of the actual time dependent coefficients (see e.g. Yu et al. (2000); Whitney (2008) for the benefit of keeping the time dependent coefficients). This leads to the Redfield Master Equation with asymptotic coefficients (Redfield Equation (asymp.)). Both variants of the Redfield Equation are not of GKSL-form.

We note in passing that for the same perturbative regime (including time dependent coefficients) a positivity preserving quantum trajectory description is available Yu et al. (1999); de Vega et al. (2005).

III.2.2 Quantum Optical Master Equation

With the aim to enforce the GKSL-form for the master equation, Eq. (8) is rewritten with $\tilde{L}(t)=\sum_{\omega}e^{\mathrm{i}\omega t}L^{\dagger}_{\omega}$ and $\tilde{L}(s)=\sum_{\omega^{\prime}}e^{-\mathrm{i}\omega^{\prime}s}L_{\omega^{\prime}}$. As before, for a sufficiently fast decaying BCF the integral can be approximated by replacing $\tilde{\rho}(s)$ with $\tilde{\rho}(t)$. The resulting equation (8) takes the form

[TABLE]

If the magnitude of $F(\omega^{\prime},t)\sim\eta/\gamma$, which represents the coupling strength, is significantly smaller than the smallest non-zero transition frequency ($\eta/\gamma\ll\min_{\omega\neq\omega^{\prime}}|\omega-\omega^{\prime}|$), it can be argued that so-called secular terms (summands with $\omega\neq\omega^{\prime}$) average to zero because of the fast oscillating phase. Keeping only the contributions $\omega=\omega^{\prime}$ and replacing $F(\omega^{\prime},t)$ by the asymptotic values $F(\omega^{\prime})=J(\omega^{\prime})+\mathrm{i}S(\omega^{\prime})$ yields, in the Schrödinger picture, the well-known Quantum Optical Master Equation of GKSL-form

[TABLE]

Note, since the so-called Lindblad Operators $L_{\omega}$ depend on the eigenvalues of $H_{\mathrm{sys}}$, for the two-qubit-system the equation changes discontinuously with the detuning of the two qubits. In the general case ($\omega_{A}\neq\omega_{B}$), the only non-zero $L_{\omega}$ are single qubit operators and read

[TABLE]

with $|{\psi_{\pm}}\rangle$ being the eigenvectors of the Pauli matrix $\sigma_{x}$ with eigenvalue $\pm 1$. In contrast, for the resonant case ($\omega_{A}=\omega_{B}=\omega$) the Lindblad operators are non local

[TABLE]

Note, for the sake of readability we write $|{\psi^{X}_{\pm}}\rangle$ instead $|{\psi_{\pm}}\rangle^{X}$. In that case additional terms appear in the master equation which are proportional to, for example, $\sigma_{z}^{A}\sigma_{z}^{B}$. Even for infinitesimally detuned qubits, these terms are missing due to the secular approximation which particularly influences the dynamics of the correlations of the two qubits (see Sec. IV.5).

III.2.3 Partial RWA

For a small detuning the unphysical discontinuity can be circumvented by using the Lindblad operators of the resonant case (Eq. (18)) also for the detuned case. Formally this corresponds to a way of deriving a master equation of GKSL kind where the RWA is applied only partially Vogt et al. (2013); Jeske et al. (2015); Tscherbul and Brumer (2015). As for the derivation of the Quantum Optical Master Equation (full RWA) Eq. (15) serves as starting point. Given that the transition frequencies $\omega$ can be grouped such that for each member $\omega$ of the group $G_{\bar{\omega}}$ the approximation $F(\bar{\omega})\approx F(\omega)$ holds, Eq. (15) becomes

[TABLE]

where $L_{\bar{\omega}}=\sum_{\omega\in G_{\bar{\omega}}}L_{\omega}$. Applying the RWA on the basis of the frequencies $\bar{\omega}$ and transforming back to the Schrödinger picture yields a master equation of GKSL form

[TABLE]

III.2.4 Coarse-Grained Master Equation

Applying a coarse-graining procedure provides yet another way to improve on the limitation of the RWA for detuned qubits while keeping the GKSL-property of the master equation Schaller and Brandes (2008); Benatti et al. (2009); Majenz et al. (2013). The method is based on a second order expansion of the time evolution operator $U(t,t+\tau)$ in the full interaction picture which yields

[TABLE]

where $\tilde{H}(s)$ is the remaining interaction Hamiltonian in the interaction picture.

Evaluating the trace over the environment on the right hand side is again done approximately by assuming that $\tilde{\rho}_{\mathrm{tot}}(t)$ can be replaced by $\tilde{\rho}(t)\tilde{\rho}_{\mathrm{env}}$ where $\mathrm{Tr}_{\mathrm{env}}[\tilde{H}(s),\tilde{\rho}(t)\tilde{\rho}_{\mathrm{env}}]=0$ has to hold333As for Eq. (7), for thermal states in combination with the usual interaction Eq. (1), this is valid.. We get

[TABLE]

This expression suggests to generate the time discrete dynamics by sequentially applying $\mathcal{Z}_{\tau}$ such that $\tilde{\rho}(t+n\tau)=\left(\mathds{1}+\mathcal{Z}_{\tau}\right)^{n}\tilde{\rho}(t)$ – provided the product state assumption is consistent at each time step. In this sense $\tau$ is related to the decay of bath correlations – on that time scale correlations between the system and the environment are expected to become unimportant for the reduced dynamics (see also Ref. Benatti et al. (2009); Majenz et al. (2013) for a discussion on the validity of the Coarse-Grained Master Equation).

However, it has been pointed out that the discrete map $\mathcal{Z}_{\tau}$ is not completely positive Schaller and Brandes (2008) – yet it is a valid GKSL generator. Therefore, if the finite difference may well be approximated by the time derivative of the reduced state, Eq. (22) turns into a master equation of GKSL-type Schaller and Brandes (2008)

[TABLE]

Note that in the mathematical limit $\tau\rightarrow 0$ the double time integral in Eq. (22) scales as $\tau^{2}$. Thus, for the Coarse-Grained Master Equation to be meaningful, a time scale separation as for the Redfield Equation and Quantum Optical Master Equation is required where the coarse-graining time has to satisfy $\tau_{\mathrm{env}}\ll\tau\ll\tau_{\mathrm{ind}}$. Again, $\tau_{\mathrm{ind}}$ is the timescale on which the reduced state changes in the interaction picture and $\tau_{\mathrm{env}}$ is the timescale set by the decay of the BCF.

To actually solve the Coarse-Grained Master Equation numerically, we do not use its formulation in obvious GKSL-form Schaller and Brandes (2008); Benatti et al. (2009). It seems more convenient to rewrite Eq. (23) solely in terms of the coupling operator decomposition $L_{\omega}$

[TABLE]

and introduce the coefficients

[TABLE]

that depend on the coarse-graining parameter $\tau$. For the Lorentzian SD given in Eq. (2) the coefficients can be evaluated explicitly. For $\omega=\omega^{\prime}$ one finds

[TABLE]

and for $\omega\neq\omega^{\prime}$

[TABLE]

As expected, when changing back to the Schrödinger picture with respect to the system, the usual Quantum Optical Master Equation is recovered Schaller and Brandes (2008) for $\tau\rightarrow\infty$

[TABLE]

We kept the Coarse-Grained Master Equation in the interaction picture in order to introduce the Lindbladian $\mathcal{Z}_{\tau}/\tau$ which can be used to construct yet another completely positive map.

III.2.5 ExpZ Map

As seen in Eq. (22), for an initial product state, the expression

[TABLE]

is correct up to second order in $t$. This motivates heuristically the completely positive ExpZ Map Majenz et al. (2013); Rivas (2017)

[TABLE]

which leads to the same short time behavior. For long times, on the other hand, $\mathcal{Z}_{\tau}$ approaches $\tau\mathcal{L}$, where $\mathcal{L}$ is the generator of the Quantum Optical Master Equation (16) in the interaction picture. Consequently, the long-time behavior of the ExpZ Map coincides with the dynamics of the Quantum Optical Master Equation. When solving the ExpZ Map as in the later examples, we directly evaluate the matrix exponential numerically for each time step.

IV Results

Our main result is the rigorous comparison of the various master equations by means of their deviation from the exact reduced dynamics. We stress that the positivity problem of the Redfield Equation is insignificant and show that the Redfield Equation (tdc) with time dependent coefficients results in the most accurate reduced dynamics. Even though the Redfield Equation (tdc) is not of GKSL-form, positivity issues of the reduced dynamics do not pose a severe problem because they show up only in a parameter regime where the approximations made are invalid. These two statements ultimately allow for the conclusion that whenever the reduced state obtained via Redfield Equation (tdc) violates positivity, the validity of any of the weak coupling approaches considered here is doubtful. Consequently, the lack of positivity preservation of the Redfield Equation (tdc) need not be seen as a shortcoming, but should rather be seen as a welcome feature. The failure to represent the true reduced dynamics cannot be detected by the positivity-preserving equations without reference to other methods.

In order to compare the various approaches, the exact dynamics (pseudo-mode method) is calculated up to a sufficiently large time $t_{\mathrm{max}}$ which depends on the coupling strength $\eta$ and the time scale of the BCF $\gamma^{-1}$ (see Fig. 1). The propagation time $t_{\mathrm{max}}$ is chosen such that the system-plus-pseudo-mode state $P(t)$ is close to the asymptotic state $P(\infty)$ for $t\geq t_{\mathrm{max}}$. More precisely, close refers to the condition for the relative difference $|P(t)-P(\infty)|/|P(\infty)|<0.01$ where the norm $|\cdot|$ denotes the Hilbert-Schmidt norm. Since we later distinguish the resonant and the detuned case with respect to the two qubit frequencies, it should be noted that the propagation time $t_{\mathrm{max}}$ obtained for the detuned case ($\omega_{A}=\Delta$ and $\omega_{B}=0.95\Delta$) is also used for the resonant case. This is justified because the relaxation towards the steady state is slower for the detuned in comparison to the resonant case.

The asymptotic state $P(\infty)$ is obtained by calculating the kernel of the Lindbladian of the truncated psudo-mode master equation. To obtain convergence with respect to the two-qubit state the truncation level of the psudo-mode is incremented by 4 until the change of the asymptotic two-qubit state is below $10^{-6}$. Therefore, the final truncation level $d$ satisfies $|\mathrm{Tr}_{\mathrm{PM}}P_{d}(\infty)-\mathrm{Tr}_{\mathrm{PM}}P_{d-4}(\infty)|/|\mathrm{Tr}_{\mathrm{PM}}P_{d}(\infty)|<10^{-6}$. The dependence of the truncation on the coupling strength $\eta$ and the BCF timescale $\gamma^{-1}$ is shown in Fig. 1. The final truncation level $d$ for the asymptotic system state is also used when propagating the pseudo-mode master equation in order to obtain the system dynamics which serves as exact reference

[TABLE]

IV.1 Error of the Master Equations

To provide error bounds independent of the initial state, we write $\rho(t)=\Lambda(t)\rho_{0}$ and use the linearity of the propagator $\Lambda(t)$. Decomposing an arbitrary initial two-qubit state into tensor products of Pauli matrices Gamel (2016) $\rho_{0}=\frac{1}{4}\sum_{\alpha,\beta}R_{\alpha,\beta}\sigma_{\alpha}\otimes\sigma_{\beta}$ with $\alpha,\beta=\{0,1,2,3\}$ and using $|R_{\alpha,\beta}|\leq 1$444$R_{\alpha,\beta}^{2}=\langle\sigma_{\alpha}\otimes\sigma_{\beta}\rangle^{2}\leq\langle\sigma_{\alpha}^{2}\otimes\sigma_{\beta}^{2}\rangle=1$, allows to bound the time dependent deviation as follows

[TABLE]

The partial deviation

[TABLE]

is calculated independently for each of the 16 combinations $\alpha,\beta$ by propagating the corresponding “initial condition” $\sigma_{\alpha}\otimes\sigma_{\beta}$ (which is a valid quantum state for $\alpha=\beta=0$ only).

To see the main features of the deviation, Fig. 2 provides an exemplary plot with three selected partial deviations and the overall sum. Three points should be noted. First, the perfect mixture as initial condition ($\epsilon_{0,0}$) yields, at the beginning, the smallest deviation, which, however, quickly reaches its asymptotic value. Second, the largest deviation occurs after a short propagation time for initial conditions related to the correlations between the two qubits ($\epsilon_{i,j}$ with $i,j=1,2,3$ corresponds to a non-zero Bloch-tensor as initial condition). And third, for the slightly detuned case, the deviation of the Redfield Equation (tdc) is by several orders of magnitude smaller as compared to the Quantum Optical Master Equation.

In order to show quantitatively how the error bound behaves while changing the coupling strength $\eta$ and correlation time $\gamma^{-1}$ we choose the maximum value of the time dependent error bound $\epsilon:=\max_{t\in[0,t_{\mathrm{max}}]}\sum_{\alpha,\beta}\epsilon_{\alpha,\beta}(t)$ as measure of accuracy. The value $\epsilon$ bounds the maximum deviation that can occur, independent of the initial state and time555This statement requires that the maximum error $\sum_{\alpha,\beta}\epsilon_{\alpha,\beta}(t)$ was reached within the time interval of propagation $[0,t_{\mathrm{max}}]$ which is ensured by choosing $t_{\mathrm{max}}$ for each combination of $(\eta,\gamma^{-1})$ such that the system has almost reached its asymptotic state (see Fig. 1).. This allows us to compare the accuracy of the various approximative methods while changing the environment.

The results are unambiguous for the detuned (Fig. 3) as well as the resonant case (Fig. 4), clearly favoring the Redfield Equation (tdc) over all other approaches considered here. Nonetheless, additional information can be drawn from these figures.

• •

The lines of constant absolute error bound can well be described by simple scaling laws in the relevant parameter regime (see Fig. 3 and Fig. 4). For the Quantum Optical Master Equation, the related ExpZ Map and the partial RWA we find an exponent $-1$, corresponding to $\eta/(\gamma\Delta)=\mathrm{const}$, as expected from a straight forward weak coupling assumption. The Redfield Equation (tdc), on the other hand, shows an exponent $-3$, corresponding to lines $\eta/(\gamma\Delta)(\Delta/\gamma)^{2}=\mathrm{const}$. The superiority of the Redfield Equation (tdc) becomes evident through the additional factor $(\Delta/\gamma)^{2}$. The lines of constant error bound for the Coarse-Grained Master Equation again follow a scaling law with exponent $-1$. Each line, however, kinks at a critical correlation time which in turn scales with the coarse graining time $\tau$ (see Fig. 8). The kink reflects an intrinsic error of the Coarse-Grained Master Equation imposed by the condition $\tau_{\mathrm{env}}\ll\tau$ which is independent of the coupling strength $\eta$.

• •

While the error bound landscape of the Quantum Optical Master Equation and the related ExpZ Map for the detuned case differs significantly from the resonant case, it hardly changes for the other methods. The explanation is found in the degeneracy of the resonant system Hamiltonian which results in the Lindblad operators $|{\psi^{A}_{-}}\rangle\langle\psi^{A}_{+}|+|{\psi^{B}_{-}}\rangle\langle\psi^{B}_{+}|$ (Eq. 18) and its Hermitian conjugate. Such Lindblad operators result in different features of the reduced dynamics as compared to the detuned case where the Lindblad operators are solely local operators of the form $|{\psi^{X}_{-}}\rangle\langle\psi^{X}_{+}|$ ($X=A,B$) and its Hermitian conjugate (see Sec. III.2.2). As a consequence of that, the general detuned-case Quantum Optical Master Equation misses some features of the dynamics of the correlations within the qubit system. More details can be found in Sec. IV.5.

• •

The error bound landscape of the ExpZ Map and the Quantum Optical Master Equation are very similar. The small advantage for the ExpZ Map can be understood by noting that the deviation of the Quantum Optical Master Equation reaches its maximum very quickly (see Fig. 2). The ExpZ Map however, yields the correct dynamics for very short times and approaches the dynamics of the Quantum Optical Master Equation for large times. Therefore the deviation of the ExpZ Map looks like the deviation of the Quantum Optical Master Equation but with a suppressed maximum at the beginning.

• •

For a wide range of environmental parameters the partial RWA master equation is the most accurate among the methods of GKSL-type. It seems to extend the levels of constant error bound of the Coarse-Grained Master Equation ($\tau\Delta=1$) beyond the kink. Notably, further increasing the detuning worsens the error of the partial RWA only little (not shown here). However, one should bear in mind that as of the particular structure of the two-qubit Hamiltonian many Lindblad operators $L_{\omega}$ vanish, which makes it obvious how to apply the partial RWA. This might not be the case for more general system Hamiltonians. As such, the partial RWA takes a special role compared to the other methods.

• •

Concerning the Coarse-Grained Master Equation ($\tau\Delta=1$) the error bound landscape is not affected by the detuning of the two qubits, just like the Redfield Equation (tdc). In contrast to the other methods, when decreasing the coupling strength only, the error bound saturates to a minimal value, which in turn depends on the correlation time. This hints again at the fact that for the Coarse-Grained Master Equation to be applicable, the correlations between the system and the environment need to become irrelevant on a faster timescale than the coarse-graining time, irrespective of the coupling strength (see also the discussion in Sec. III.2.4). However, for the detuned case Fig. 3 shows that there is a regime (small correlation time and fairly large coupling strength) where the Coarse-Grained Master Equation is more accurate than the Quantum Optical Master Equation and the ExpZ Map (for more details see Sec. IV.4 and Sec. IV.5).

The discussion so far has ignored the main criticism concerning the Redfield Equation (tdc), the lack of guaranteed positivity. By choosing a physical state as initial condition ($\psi_{0}=|{\uparrow\uparrow}\rangle$) we are able to keep track of the positivity of the reduced state. Further, for a particular initial state the relative error $r(t)=|\rho_{\mathrm{ref}}(t)-\rho_{\mathrm{M}}(t)|/|\rho_{\mathrm{ref}}(t)|$ can be calculated which allows a comparison of the methods based on the actual error instead of the error bound used earlier. Nonetheless, since it turns out that the relative error landscape for each method is very similar to the error bounds shown in Fig. 3 and Fig. 4 (therefore it is not shown here) the initial condition $\psi_{0}=|{\uparrow\uparrow}\rangle$ can be seen as a generic initial condition, not featuring any special behavior with respect to the applicability of the various master equations.

In Fig. 5 the parameter region where the maximum relative error is below $5\%$ is shown for the more challenging case of two detuned qubits ($\omega_{A}=\Delta$ and $\omega_{B}=0.95\Delta$) with initial condition $\psi_{0}=|{\uparrow\uparrow}\rangle$. The earlier picture from the initial state independent discussion is restated: the Redfield Equation (tdc) covers the largest parameter region followed by the Redfield Equation (asymp.). The ExpZ Map performs slightly better than the Quantum Optical Master Equation. The Coarse-Grained Master Equation with $\tau\Delta=1$ outreaches the ExpZ Map for sufficiently short correlations times, however, is outperformed by the partial RWA. To add to this picture, keeping track of the positivity for the reduced dynamics obtained from the Redfield Equation (tdc) reveals that positivity problems do only occur in a parameter region where the Redfield Equation (tdc) becomes significantly invalid. One can even go further by reading the plot in Fig. 5 such that a positivity violation of the reduced dynamics obtained from the Redfield Equation (tdc) allows to keep track of the validity of the underlying approximations made, without having to refer to the exact solution.

Thus, the criticism directed at the Redfield Equation for not being of GKSL-form may be refuted considerably in the light of its accuracy and, in particular, the benefit of using the positivity violation of the Redfield Equation (tdc) as a criterion for any of the perturbative master equations considered here to be applicable.

IV.2 The Advantage of Time Dependent Coefficients

It should be emphasized that the error of the Redfield Equation (asymp.)  with asymptotic coefficients is slightly larger than the error of the Redfield Equation (tdc)  with time dependent coefficients (see Fig. 6). However, the Redfield Equation (asymp.) still outperforms the other methods under consideration. Notably, even in a regime where the relative error is fairly small, transiently non-positive reduced states may occur when using the asymptotic coefficients. Of course, the order of magnitude of the negative eigenvalue does not exceed the order of the error (see Fig. 6).

The difference between the two variants of the Redfield Equation is shown in Fig. 6, where the maximum relative difference and the minimum negative eigenvalue of the dynamics are plotted for a slice through the parameter space with fixed coupling strength. Although only small in magnitude, non-positive eigenvalues of the Redfield Equation (asymp.) dynamics occur already for correlation times where the relative error is still small. When increasing the correlation time, the non-positive eigenvalues increase in magnitude roughly in the same manner as the relative error. In contrast, for the Redfield Equation (tdc) the non-positivity sets in suddenly.

Examining the time dependence of the smallest eigenvalue (see the right panels in Fig. 6) suggests that there are two causes for the positivity violation. First, using asymptotic coefficients as in the Redfield Equation (asymp.), obviously, is not justified for the initial dynamics on the time scale of the correlation time. As a result, non-positive eigenvalues occur during that initial dynamics. Their magnitude decreases with decreasing correlation time which is in line with the observation that for a delta-like correlation function using the asymptotic coefficients becomes exact. However, the non-positive eigenvalues occurring during the initial dynamics disappear after the correlation time has passed (this initial positivity problem is often discussed in terms of an initial slippage Haake and Lewenstein (1983); Suárez et al. (1992); Gaspard and Nagaoka (1999); Yu et al. (2000); Cheng and Silbey (2005)). Using the time dependent coefficients as in the Redfield Equation (tdc) circumvents this problem entirely (see the useful Ref. Whitney (2008) for a thorough investigation of this phenomenon with analytical results for very short correlation times).

The second reason simply originates from the fact that for larger correlation times (or larger coupling strengths) the perturbative approach of the Redfield Equation in general (both time dependent coefficients and asymptotic rates) becomes invalid, resulting in long lasting violation of the positivity (and accuracy) of the reduced dynamics.

IV.3 Linear Scaling of the Error

Concerning the scaling of the error with the coupling strength it has been shown that a perturbative master equation of order $2n$ in the coupling strength yields an accuracy for the long time dynamics which is of the order $2n-2$ Fleming and Cummings (2011). Therefore for the second order master equations considered here the scaling has to be as good as zeroth order and cannot be, in general, of second order. On the other hand, since it is also known that the Quantum Optical Master Equation becomes exact in the zero coupling (scaling) limit Davies (1974) the error has to vanish (at least for the Quantum Optical Master Equation). In Fig. 7 the scaling of the error is shown for different environmental correlation times $\Delta/\gamma$. For all of them the plots suggest a linear behavior for the Redfield Equation (tdc), Redfield Equation (asymp.), Quantum Optical Master Equation and ExpZ Map. However, in the case of the Coarse-Grained Master Equation and the partial RWA the error seems to decrease as well until it reaches a finite value. In case of the partial RWA this method intrinsic error originates from replacing $F(\omega)$ with the corresponding value $F(\bar{\omega})$ for the cluster frequency $\bar{\omega}$. This remaining error is consistent, since in the zero coupling limit the Quantum Optical Master Equation does distinguish even between very close transitions frequencies. For the Coarse-Grained Master Equation the condition $\tau_{\mathrm{env}}\ll\tau$ has to be met which also induces a coupling strength independent contribution to the overall error (see also next Sec. IV.4).

IV.4 The Coarse-Graining Time $\tau$

From a mathematical point of view, the coarse-graining parameter $\tau$ can be chosen freely. However, we have already stressed in Sec. III.2.4 (see also Ref. Benatti et al. (2009); Majenz et al. (2013)) that in order to relate the resulting dynamics to the microscopic model, $\tau$ has to fulfill two conditions. By physical means the condition $\tau_{\mathrm{env}}\ll\tau$ justifies the product state replacement of the total state $\tilde{\rho}(\tau)\rightarrow\tilde{\rho}_{\mathrm{sys}}(\tau)\otimes\tilde{\rho}_{\mathrm{env}}$ after the first time step $\tau$ and, thus, allows to iteratively propagate subsequent time steps $\tau$ Majenz et al. (2013). The other condition $\tau\ll\tau_{\mathrm{ind}}$, where $\tau_{\mathrm{ind}}\sim\gamma/\eta$, ensures sufficiently slow system dynamics in the interaction picture, such that the finite difference is well represented by the derivative Benatti et al. (2009); Majenz et al. (2013).

Notably, the time scale set by the energy differences of the system Hamiltonian does not play a role. Consequently, for suitable environments, where the above time scale separation holds, the Coarse-Grained Master Equation is applicable irrespectively of the system Hamiltonian and, thus, provides a master equation beyond the RWA.

The influence of the coarse-graining parameter $\tau$ on the error landscape is shown in Fig. 8. To examine the effect of the “coarse-graining” of the Coarse-Grained Master Equation we also show the error landscape where the exact reduced state averaged over the coarse-graining parameter $\tau$ in the interaction picture

[TABLE]

is used as reference. The average ensures that $\langle\tilde{\rho}_{\mathrm{ref}}(0)\rangle_{\tau}=\rho_{\mathrm{ref}}(0)$ which also serves as initial condition for the Coarse-Grained Master Equation. Fig. 8 shows that there is no significant difference in the overall error behavior between the two cases (upper row: $\tilde{\rho}_{\mathrm{ref}}(t)$, lower row $\langle\tilde{\rho}_{\mathrm{ref}}(t)\rangle_{\tau}$). However, minor differences can be noted in the regime where the error is already small. In that case, the dynamics obtained from the Coarse-Grained Master Equation matches the $\tau$-averaged exact dynamics better than the non-averaged exact dynamics.

Additionally, the plots in Fig. 8 show explicitly that for a very small correlation time, such that a rather small coarse-graining time is justified, the Coarse-Grained Master Equation is also applicable for somewhat stronger couplings, a regime in general not accessible by the Quantum Optical Master Equation. This statement will become more explicit in the example dynamics shown in the following.

IV.5 Influence of the Secular Approximation on the Qubit Correlations

Recall, for the detuned case the Lindblad operators read $L_{\omega_{A}}=\frac{1}{2}|{\psi_{-}}\rangle^{A}\langle\psi_{+}|^{A}$, $L_{\omega_{B}}=\frac{1}{2}|{\psi_{-}}\rangle^{B}\langle\psi_{+}|^{B}$ and their Hermitian conjugate. Viewing $\omega_{B}$ as a parameter of the corresponding Quantum Optical Master Equation (fix the form of the Lindblad operators), the resonant case $\omega_{A}=\omega_{B}$ can also be treated with that Quantum Optical Master Equation. On the other hand, for the resonant case the Lindblad operators $L_{\omega}=\frac{1}{2}\left(|{\psi_{-}}\rangle^{A}\langle\psi_{+}|^{A}+|{\psi_{-}}\rangle^{B}\langle\psi_{+}|^{B}\right)$ and its Hermitian conjugate can be derived explicitly, resulting in a different master equation of GKSL-form. The difference of the two variants becomes obvious by realizing that the Lindblad operators enter the master equation quadratically. For example, the Lamb-shift Hamiltonian for the Lindblad operators derived from the detuned case, however used in resonance $\omega=\omega_{A}=\omega_{B}$, reads

[TABLE]

In contrast, using the Lindblad operators $L_{\omega}=\frac{1}{2}\left(|{\psi^{A}_{-}}\rangle\langle\psi^{A}_{+}|+|{\psi^{B}_{-}}\rangle\langle\psi^{B}_{+}|\right)$ additional terms occur in the Lamb-shift Hamiltonian

[TABLE]

These additional terms, which effectively result in a unitary coupling between the two qubits Benatti et al. (2003); Tana and Ficek (2004), are missing due to the RWA applied in the detuned case. In the same manner, differences between the two variants of the Quantum Optical Master Equation occur also in the dissipator. The non-local structure (in terms of the two qubits) of the additional contribution will particularly influence the dynamics of the two-qubit correlations (see Fig. 9 and Fig. 10).

To summarize, the special Quantum Optical Master Equation derived for the resonance condition includes non-local terms expected to influence the correlation dynamics of the qubits. Once the detuned case is considered, the formalism of the Quantum Optical Master Equation results in an equation without such non-local terms. It is precisely the motivation of the partial RWA, Coarse-Grained Master Equation and the ExpZ Map to overcome this shortcoming Schaller and Brandes (2008); Benatti et al. (2009, 2010); Majenz et al. (2013); Rivas (2017).

In order to exemplify how the various approaches approximate the dynamics, we pick two pairs of $\eta/\Delta^{2}$ and $\Delta/\gamma$ where the differences are sufficiently well visible. In particular we distinguish between the dynamics of the local expectation value $\langle\mathds{1}\otimes\sigma_{z}\rangle$ and the non-local quantity $\langle\sigma_{z}\otimes\sigma_{z}\rangle$.

In Fig. 9 the dynamics for a rather weak coupling strength $\eta/\Delta^{2}=0.149$ and a correlation time $\Delta/\gamma=0.673$, which is of the order of the single qubit time scale, is shown. For the resonant case, all methods except the Coarse-Grained Master Equation approximate the exact dynamics very well. As expected, for the slightly detuned case, where the detuning results in an additional system time scale slower than the correlation time, the validity of the Quantum Optical Master Equation breaks down. However, the single qubit dynamics is well recovered. Significant deviations are visible for the correlation dynamics of the two qubits. The ExpZ Map smoothly interpolates from the exact dynamics for short times to the values of the Quantum Optical Master Equation for longer times. Concerning the Coarse-Grained Master Equation, the difference to the exact dynamics is equally visible for both, the local and non-local expectation value independently of the detuning. This is plausible, because the coarse-graining time $\Delta\tau=1$ is of the same order as the correlation time $\Delta/\gamma=0.673$ which renders the Coarse-Grained Master Equation to be inaccurate. The Redfield Equation (tdc) and the partial RWA, however, can hardly be distinguished and match the exact dynamics even for the non-local contribution in the detuned case.

In the next example, the coupling strength is chosen larger $\eta/\Delta^{2}=1.29$ while the correlation time becomes shorter $\Delta/\gamma=0.165$. Again, as expected, for the resonant case all methods yield suitable dynamics for the local as well as the non-local expectation value. In the detuned case this holds for the local quantity, too. However, both, the Quantum Optical Master Equation and the ExpZ Map, do not account for the slow decay of the $\sigma_{z}\otimes\sigma_{z}$ correlations. In contrast, as of the shorter correlation time, the Coarse-Grained Master Equation is more suitable as in the previous example and, thus, particularly outperforms the Quantum Optical Master Equation and ExpZ Map on the correlation dynamics. The same holds true for the partial RWA. Again, for all examples, the Redfield Equation (tdc) provides the most accurate results.

In all examples discussed so far the partial RWA method yields more accurate results than the Coarse-Grained Master Equation. However, as shown in the previous Sec. IV.4, for sufficiently small correlations times $\Delta/\gamma$ a coarse-graining time $\tau\Delta<1$ is justified which is beneficial in particular for large coupling strength. In Fig. 11 we show that in this regime the Coarse-Grained Master Equation with $\tau\Delta=0.3$ is significantly more accurate than the partial RWA.

V Conclusions

The need to describe the dynamics of open quantum systems most adequately has led to a wealth of perturbative master equations. While each approach was developed with a certain objective in mind a general and comparative assessment of the accuracy of all these methods was missing. With our work we fill this gap by quantitatively comparing the approximate reduced state with the exact dynamics. Notably, whether or not a particular master equation guarantees positivity we consider secondary. Our error-based examination confirms that whenever a perturbative approach is justified the Redfield Equation (tdc) is the method of choice. As indicated by an initial state independent error bound, the Redfield Equation (tdc) substantially outperforms the other methods considered here (Redfield Equation (asymp.), Quantum Optical Master Equation, ExpZ Map, Coarse-Grained Master Equation and partial RWA). Further, the lack of ensured positivity preservation should not be considered as a bug, but as a feature: it indicates the breakdown of the weak coupling approximation.

In order to contribute to a better understanding of the applicability of the various master equations, we have also investigated their error in detail. For the Quantum Optical Master Equation we have explicitly argued – and confirmed by examples – that in the general detuned case of the two qubits, the RWA most significantly effects the correlations between the two qubits. The approaches ExpZ Map, partial RWA and Coarse-Grained Master Equation do – to some extent – improve on the shortcomings of the Quantum Optical Master Equation as they do not explicitly make use of the RWA while yielding positive dynamics. Our error analysis reveals that the ExpZ Map performs slightly better in terms of the maximum error for the entire dynamics but mimics qualitatively the same error landscape as the Quantum Optical Master Equation. Moreover, we find that whenever the time scale separation $\tau_{\mathrm{env}}\ll\tau\ll\tau_{\mathrm{ind}}$ is satisfied, the Coarse-Grained Master Equation yields good results, irrespectively of the system Hamiltonian, that is, it does not distinguish between the detuned and resonant case. The coarse-graining parameter $\tau$-dependent error landscape also qualitatively differs from the Quantum Optical Master Equation correspondent. Exploiting this feature allowed us to explicitly show that there is a region in the parameter space spanned by the coupling strength $\eta/\Delta^{2}$ and correlation time $\Delta/\gamma$ where the Coarse-Grained Master Equation outperforms the Quantum Optical Master Equation significantly. The partial RWA takes a special role since the applicability depends on specific spectral features of the system Hamiltonian. For the non-interacting two qubit system, considered here, with only two relevant frequencies the implementation is straight forward. For that case our analysis reveals that in most cases the partial RWA is more accurate than the other methods of GKSL-kind. Only when smaller $\tau$ are justified the partial RWA is outperformed by the Coarse-Grained Master Equation in the short correlation time and large coupling strength regime.

Although we focus on a particular system of two qubits and a Lorentzian environment, we are confident that our conclusions hold true for generic systems that contain a wide range of transition frequencies.

VI Acknowledgments

Fruitful discussions with Kimmo Luoma and Sebastian Diehl are gratefully acknowledged. The computations were performed on a Bull Cluster provided at the Center for Information Services and High Performance Computing (ZIH) at TU Dresden. Support by the IMPRS at the Max Planck Institute for the Physics of Complex Systems (Dresden) and in part by the National Science Foundation under Grant No. NSF PHY-1748958 during time at KITP (UCSB) is gratefully acknowledged. We thank Elke and Steve Langdon: “There is too much negativity surrounding the issue of positivity”.