Correlation Picture Approach to Open-Quantum-System Dynamics

TL;DR

This paper introduces a correlation picture for open quantum systems, deriving an exact Lindblad-like dynamical equation that includes initial correlations, and develops a weak-correlation expansion for systematic approximations, improving over existing methods.

Contribution

The paper presents a novel correlation picture framework that yields an exact open-system dynamics equation including initial correlations, along with a systematic perturbative expansion for practical approximations.

Findings

Exact Lindblad-like equation with initial correlations derived

Weak-correlation expansion formalism developed for approximations

New Markovian master equation differing from existing approaches

Abstract

We introduce a new dynamical picture, referred to as correlation picture,' which connects a correlated state to its uncorrelated counterpart. Using this picture allows us to derive an exact dynamical equation for a general open-system dynamics with system--environment correlations included. This exact dynamics is in the form of a Lindblad-like equation even in the presence of initial system-environment correlations. For explicit calculations, we also develop a weak-correlation expansion formalism that allows us to perform systematic perturbative approximations. This expansion provides approximate master equations which can feature advantages over existing weak-coupling techniques. As a special case, we derive a Markovian master equation, which is different from existing approaches. We compare our equations with corresponding standard weak-coupling equations by two examples, where our correlation picture formalism is more accurate, or at least as accurate as weak-coupling equations.

Figures

Equations

$${\varrho }_{\mathsf{SB}}(\tau )={\varrho }_{\mathsf{S}}(\tau )\otimes {\varrho }_{\mathsf{B}}(\tau )+\chi (\tau ),$$

$$\chi =-i[H_{\chi },{\varrho }_{\mathsf{S}}\otimes {\varrho }_{\mathsf{B}}].$$

$$\chi =-i[[H_{\chi },{\varrho }_{\mathsf{S}}\otimes {\varrho }_{\mathsf{B}}]].$$

$${\varrho }_{\mathsf{SB}}=\text{mathpzc}{E}_{\chi }[{\varrho }_{\mathsf{S}}\otimes {\varrho }_{\mathsf{B}}]:={\varrho }_{\mathsf{S}}\otimes {\varrho }_{\mathsf{B}}-i[[H_{\chi },{\varrho }_{\mathsf{S}}\otimes {\varrho }_{\mathsf{B}}]].$$

$$\text{mathpzc}{D}_{\text{mathpzc}s}[{\varrho }_{\mathsf{SB}}(\tau )]=$$

$$\displaystyle\mathpzc{D}_{\mathpzc{c}}[\circ]=-i[H_{\mathsf{SB}},\circ]-\big{[}H_{\mathsf{SB}},[\hskip-2.56073pt[H_{\chi},\circ]\hskip-2.56073pt]\big{]}.$$

$${\dot{\varrho }}_{\mathsf{S}}=$$

$$\displaystyle-\textstyle{\sum_{i\neq 0,j}}\big{[}\mathpzc{S}_{i},c_{ij}\mathpzc{S}_{j}\varrho_{\mathsf{S}}-c_{ij}^{\ast}\varrho_{\mathsf{S}}\mathpzc{S}_{j}\big{]},$$

$$c_{ij}(\tau )=⟨\text{mathpzc}{B}_i\text{mathpzc}{B}_j^{\chi }(\tau ){⟩}_{\mathsf{B}}$$

$${\dot{\varrho }}_{\mathsf{S}}$$

$$\displaystyle=-i[H_{\mathsf{S}}+{\mathbbmss{h}}_{{\mathsf{L}}}^{\chi},\varrho_{\mathsf{S}}]+\sum_{m}\gamma_{m}^{\chi}\big{(}2L^{\chi}_{m}\varrho_{\mathsf{S}}L^{\chi{\dagger}}_{m}-\{L^{\chi{\dagger}}_{m}L^{\chi}_{m},\varrho_{\mathsf{S}}\}),$$

$${\dot{\varrho }}_{\mathsf{S}}=$$

$$-{\gamma }_2^{\chi }(2{\sigma }_{+}{\varrho }_{\mathsf{S}}{\sigma }_{-}-\{{\sigma }_{-}{\sigma }_{+},{\varrho }_{\mathsf{S}}\}),$$

$$c_{ij}(\tau )$$

$$a_{ij}(\tau )\approx$$

$${\text{mathbbmss}h}_{\mathsf{L}}^{\chi }(\tau )\approx$$

$$-2\tau {\sum }_{(i,j)\ne (0,0)}⟨\text{mathpzc}{S}_j{⟩}_{{\mathsf{S}}_0}\mathrm{Im}⟨\text{mathpzc}{B}_i\text{mathpzc}{B}_j{⟩}_{{\mathsf{B}}_0}\text{mathpzc}{S}_i,$$

$$\mathbf{\chi }(\tau )={\sum }_{\ell =0}^{\infty }{\mathbf{f}}_{\ell }(\tau ),$$

$${\varrho }_{\mathsf{SB}}(\tau +\xi )=U_{\xi }{\varrho }_{\mathsf{S}}(\tau )\otimes {\varrho }_{\mathsf{B}}(\tau )U_{\xi }^{†}+U_{\xi }\chi (\tau )U_{\xi }^{†},$$

$$\dot{\mathbf{\chi }}(\tau )=-i[{\tilde{\mathbf{H}}}_{\mathrm{I}}(\tau ),{\mathbf{\varrho }}_{\mathsf{S}}(\tau )\otimes {\mathbf{\varrho }}_{\mathsf{B}}(\tau )]+{\mathcal{Y}}_{\tau }[\mathbf{\chi }(\tau )],$$

$${\mathcal{Y}}_{\tau }[\circ ]:=$$

$$+i{\mathbf{\varrho }}_{\mathsf{S}}(\tau )\otimes {\mathrm{Tr}}_{\mathsf{S}}[{\mathbf{H}}_{\mathrm{I}}(\tau ),\circ ].$$

$$\mathbf{\chi }(\tau )=$$

$$\displaystyle+\textstyle{\sum_{k=0}^{\infty}}\int_{0}^{\tau}ds_{1}\int_{0}^{s_{1}}ds_{2}\,\cdots\int_{0}^{s_{k-1}}ds_{k}\,\bm{\mathcal{Y}}_{s_{1}}\left[\bm{\mathcal{Y}}_{s_{2}}\cdots\Big{[}\bm{\mathcal{Y}}_{s_{k}}\big{[}-i\int_{0}^{s_{k}}ds\,[\widetilde{\bm{H}}_{\mathrm{I}}(s),\bm{\varrho}_{\mathsf{S}}(s)\otimes\bm{\varrho}_{\mathsf{B}}(s)]\big{]}\Big{]}\ldots\right],$$

$${\mathbf{\chi }}^{(1)}(\tau )=-i{\int }_0^{\tau }ds[{\tilde{\mathbf{H}}}_{\mathrm{I}}(s),{\mathbf{\varrho }}_{\mathsf{S}}(s)\otimes {\mathbf{\varrho }}_{\mathsf{B}}(s)].$$

$${\mathbf{\chi }}_{\mathrm{MLL}}^{(1)}(\tau )\approx -i\tau [{\tilde{\mathbf{H}}}_{\mathrm{I}}(\tau ),{\mathbf{\varrho }}_{\mathsf{S}}(\tau )\otimes {\mathbf{\varrho }}_{\mathsf{B}}(\tau )],$$

$$H_{\mathsf{SB}}={\omega }_0{\sigma }_{+}{\sigma }_{-}+{\sum }_n{\omega }_n{\hat{a}}_n^{†}{\hat{a}}_n-{\sigma }_x\otimes \text{mathpzc}{O}_{\mathsf{B}},$$

$$\dot{\varrho}_{\mathsf{S}}(\tau)=-i[H_{\mathsf{S}},\varrho_{\mathsf{S}}(\tau)]+\gamma(\tau)\big{(}\sigma_{x}\varrho_{\mathsf{S}}(\tau)\sigma_{x}-\varrho_{\mathsf{S}}(\tau)\big{)},$$

$$\displaystyle\varrho_{\mathrm{ee}}(\tau)=1/2+\big{(}\varrho_{\mathrm{ee}}(0)-1/2\big{)}e^{-2\tau^{2}\,\mathrm{Cov}_{\mathsf{B}_{0}}(\mathpzc{O}_{\mathsf{B}},\mathpzc{O}_{\mathsf{B}})}.$$

$$H_{\mathsf{SB}}={\omega }_0{\hat{a}}^{†}\hat{a}+{\sum }_{k=1}^M{\omega }_k{\hat{b}}_k^{†}{\hat{b}}_k+{\sum }_{k=1}^M{g}_k({\hat{a}}^{†}{\hat{b}}_k+\hat{a}{\hat{b}}_k^{†}),$$

$${\dot{\varrho }}_{\mathsf{S}}(\tau )=$$

$$\displaystyle-\varrho_{\mathsf{S}}(\tau)\hat{a}^{{\dagger}}\hat{a}\big{)}.$$

$${\dot{\varrho }}_{\mathsf{S}}(\tau )=$$

$$\displaystyle+\gamma\big{(}2\hat{a}\varrho_{\mathsf{S}}(\tau)\hat{a}^{{\dagger}}-\hat{a}^{{\dagger}}\hat{a}\varrho_{\mathsf{S}}(\tau)-\varrho_{\mathsf{S}}(\tau)\hat{a}^{{\dagger}}\hat{a}\big{)},$$

$$\varrho^{{\dagger}}_{\mathsf{S}}(\tau)\otimes\varrho^{{\dagger}}_{\mathsf{B}}(\tau)\big{[}iH^{{\dagger}}_{\chi}(\tau)\big{]}+\big{[}iH^{{\dagger}}_{\chi}(\tau)\big{]}^{{\dagger}}\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)=\chi(\tau),$$

Full text

I introduction

The recent rise in high-fidelity quantum technological devices has necessitated detailed understanding of open quantum systems and how their environment influences their performance. The theory of open quantum systems has spurred numerous possibilities to harness the power of the environments in various physical tasks eng-dis ; light-harvesting . Although in this framework quantum correlations play a key role, it has remained unresolved how to explicitly keep track of correlations between a system and its environment in the dynamical equation.

The description of the joint evolution of a quantum system and its environment (or bath) through the unitary evolution given by the Schrödinger equation is hampered by the large dimensionality of the Hilbert space of the bath. Although there already exist a plethora of approximate methods book:Breuer-Petruccione ; book:Rivas-Huelga ; Nakajima ; Zwanzig ; Breuer-proj ; devega ; Maniscalco ; Breuer-Kappler for obtaining dynamics of an open quantum system by a closed set of equations, these approaches in general do not incorporate correlations between the system and its environment (see, e.g., Refs. Nazir-etal ; pollock-1 ; paz-silva for attempts to take correlations into account). Despite previous attempts to prove the existence of universally valid time-local Lindblad-like master equations for general dynamics Lindblad ; GKS ; Andersson ; Lidar-Whaley , a microscopic derivation which incorporates correlations within the dynamics has been missing. Whether a time-local master equation exists that can generally describe dynamics of systems with correlations is still an open question PRL-Chruscinski .

Here we resolve this issue by introducing the correlation picture, a new technique through which we relate any correlated state of a composite system in the Schrödinger picture to an uncorrelated description of that system. Using the dynamical picture allows us to assign a correlation parent operator to the instantaneous system-bath correlation. This correlation parent operator is used as a key element in a microscopic derivation of the dynamical equation of the system. This introduced framework for studying open quantum dynamics enables us to derive a universal Lindblad-like (ULL) equation which is time-local and valid for any quantum dynamics. The ULL form is derived without approximations and it is valid also when the system is initially correlated with the bath. A Lindblad-like form for general dynamical equations has also been noted in recent attempts to generalize quantum mechanics Weinberg-Lindblad ; Bassi-Lindblad ; Diosi-Lindblad .

We study the behavior of the ULL equation under chosen approximations and are able to derive conveniently solvable master equations which to a good degree of accurately reproduce the exact dynamics in the corresponding parameter regimes. In particular, in the vicinity of time instants where the correlations become negligible, the ULL equation reduces to a Markovian Lindblad-like (MLL) master equation, in which the jump rates are positive. We prove that this equation correctly characterizes the universal quadratic short-time behavior of the system dynamics Braun , in contrast to the standard Lindblad equation which gives a purely linear behavior for short times Rivas-short-time ; DelCampo . In addition, we demonstrate that our MLL equation, which does not utilize the secular approximation, may in some cases faithfully describe the long-time behavior of the system. This MLL equation thus constitutes a useful framework for studying open-quantum-system dynamics beyond the weak-coupling regime. Based on the MLL methodology, for explicit or practical calculations we also develop a systematic perturbative weak-correlation expansion of higher orders which include non-Markovian effects. Interestingly, we show that even the lowest order of the constructed ULL equation (ULL$2$) can outperform corresponding weak-coupling master equations. This shows that giving the principal role to correlations rather than coupling offers a strong alternative to existing weak-coupling techniques.

II Correlating transformation

Any given system-bath state at an arbitrary instant of time $\varrho_{\mathsf{SB}}(\tau)$ can be decomposed in terms of an uncorrelated part, given by the tensor product of the instantaneous reduced states of the subsystems $\varrho_{\mathsf{S}}(\tau)=\mathrm{Tr}_{\mathsf{B}}[\varrho_{\mathsf{SB}}(\tau)]$ and $\varrho_{\mathsf{B}}(\tau)=\mathrm{Tr}_{\mathsf{S}}[\varrho_{\mathsf{SB}}(\tau)]$, and the remainder $\chi(\tau)$ which carries all correlations within the total state,

[TABLE]

where $\mathrm{Tr}_{\mathsf{S}}[\chi]$ and $\mathrm{Tr}_{\mathsf{B}}[\chi]$ are null operators on the bath and system Hilbert spaces, respectively. Below, we refer to $\chi(\tau)$ as the correlation operator or simply the correlation. It includes all kinds of correlations, whether classical or quantum mechanical. The latter, in the form of entanglement, discord, or other more complex types, have a rich and resourceful nature in physics, e.g., in energy fluctuations of thermodynamical systems Sampaio and in quantum information tasks book:Nielsen ; vedral-rmp ; modi .

To define our correlation picture, we introduce an operation $\mathpzc{E}_{\chi}$ which transforms the uncorrelated state $\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)$ to the correlated state $\varrho_{\mathsf{SB}}(\tau)$. We refer to the opposite operation relating the correlated state to the uncorrelated one as decorrelating. These are interesting operations, which also appear in the context of quantum statistical physics, where they are dubbed the quantum Boltzmann map and relate to the stosszahlansatz Wolf-lectures . Decorrelating transformations have already been investigated in the literature D'Ariano-Demkowicz:PRL , and it is known that a universal decorrelating machine would violate linearity of quantum mechanics Terno . Our correlating transformation, indeed, is not universal, i.e., $\mathpzc{E}_{\chi}$ depends on the states. To gain insight on how to find $\mathpzc{E}_{\chi}$, it is helpful to start with the case of creating a weakly correlated state. We emphasize that this case study will serve merely to set the scene, whereas later we do not assume any weak correlation in our general analysis.

*Correlating transformation for a weakly correlated state.—*Let us aim to transform an uncorrelated state $\varrho_{\mathsf{SB}}=\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$ to a weakly correlated state $\varrho_{\mathsf{SB}}=\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}+\chi$ with $\|\chi\|\ll 1$ ($\|\cdot\|$ being the operator norm). A natural way to do so is to apply an entangling gate on the uncorrelated state. Consider, e.g., an entangling gate entangler3 described by a unitary transformation $U_{\chi}=e^{-iH_{\chi}}$, with $\|H_{\chi}\|\ll 1$ nb . That is, here the transformation ($U_{\chi}$) is assumed to be given and we look for the corresponding generated correlation operator $\chi$. This gate results in a weakly correlated state $U_{\chi}\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}U_{\chi}^{{\dagger}}\approx\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}-i[H_{\chi},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]$. Comparing this state with $\varrho_{\mathsf{SB}}=\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}+\chi$, we observe that up to $O(\|H_{\chi}\|^{2})$ the correlation $\chi$ obeys the following equation:

[TABLE]

We refer to the dimensionless operator $H_{\chi}$ as the correlation parent operator. The equivalent correlating transformation in this regime is obtained by inserting Eq. (2) into Eq. (1) as $\varrho_{\mathsf{SB}}=\mathpzc{E}_{\chi}[\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]:=\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}-i[H_{\chi},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]$.

Correlating transformation for a general state.—In the above scenario, we applied a known unitary transformation to create a correlated state and obtained the correlation $\chi$ in terms of the input uncorrelated state $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$. Now, we return to our question of interest: to find the operation which transforms a given $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$ to its associated $\varrho_{\mathsf{SB}}$ with a definite correlation operator $\chi$ which is not necessarily small.

Although for any given pair of quantum states $\varrho_{1}$ and $\varrho_{2}$, it is possible to find a quantum map or channel $\mathpzc{E}$ such that $\mathpzc{E}[\varrho_{1}]=\varrho_{2}$ Rabitz , such a map does not necessarily have a unitary representation. Hence if we postulate the form given in Eq. (2) as the equation connecting $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$ and $\chi$, for some $H_{\chi}$, it does not generally have a Hermitian solution for $H_{\chi}$. However, with the insight gained from the weakly correlated case, we still choose $\chi$ as in Eq. (2) but relax the Hermiticity constraint on $H_{\chi}$. Since the left side of Eq. (2) is Hermitian, to ensure the Hermiticity of the right side, with a non-Hermitian $H_{\chi}$, we introduce a generalized commutator $[\hskip-2.56073pt[A,B]\hskip-2.56073pt]=AB-B^{{\dagger}}A^{{\dagger}}$ MR and replace Eq. (2) with

[TABLE]

We consider this relation as the general definition of $H_{\chi}$, valid for any $\chi$. Fortunately, this equation has always a solution for $H_{\chi}$ provided that $P_{0}(\tau)\chi(\tau)P_{0}(\tau)=0$, where $P_{0}(\tau)$ is the projector onto the null-space of $\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)$ Djordjevic . We prove that this condition is always satisfied SM —and we also provide the solution for $H_{\chi}$. Using Eq. (3), we define our correlating transformation $\mathpzc{E}_{\chi}$ as

[TABLE]

We note that neither the solution of Eq. (3) for $H_{\chi}$ nor the formal connection of $\chi$ to $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$ is unique (one may suggest other forms). This means that there is a kind of gauge freedom in choosing the correlating transformation. However, this nonuniqueness does not affect our formalism, because any solution of Eq. (3) and any properly chosen form for $\mathpzc{E}_{\chi}$ must generate the same given $\chi$—which is the physical quantity which should remain unchanged under these gauge transformations.

Note that the uncorrelated state $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$ is not the state of the total system (because in general $\chi\neq 0$); rather, we take this state as the description of the total system in a pertaining correlation picture. In order to keep the dynamics of the state in this picture faithful to the Schrödinger equation, we need to devise an appropriate formulation. Figure 1 depicts the correlating transformation and the correspondingly emerging correlation picture—which is explained in detail below.

III Correlation picture dynamics

We aim to apply our correlation picture transformation between the correlated and uncorrelated states, $\varrho_{\mathsf{SB}}$ and $\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}$, respectively, to obtain a dynamical equation for $\varrho_{\mathsf{S}}$. Our approach can be considered in the spirit of the derivation of the Nakajima–Zwanzig (NZ) equation Nakajima ; Zwanzig . However, rather than applying a decorrelating projector to omit system–bath correlations (while implicitly carrying correlations into another complementary equation), we employ our correlating transformation within the Schrödinger equation of the total system. Hence we shall explicitly retain contributions to the system dynamics from the correlations in the total system.

Let us assume that the total Hamiltonian of the system and the bath is given by $H_{\mathsf{SB}}=H_{\mathsf{S}}+H_{\mathsf{B}}+H_{\mathrm{I}}$, where the last term denotes the system–bath interaction. We employ the correlating transformation (4) to obtain a counterpart for the Schrödinger picture generator $\mathpzc{D}_{\mathpzc{s}}[\circ]=-i[H_{\mathsf{SB}},\circ]$ (throughout this paper we have assumed the natural units $\hbar\equiv k_{\mathrm{B}}\equiv 1$). More precisely, we define this operator $\mathpzc{D}_{\mathpzc{c}}$ such that

[TABLE]

By inserting the correlating transformation (4) in the Schrödinger equation as \mathpzc{D}_{\mathpzc{s}}[\varrho_{\mathsf{SB}}(\tau)]=\mathpzc{D}_{\mathpzc{s}}\big{[}\mathpzc{E}_{\chi}[\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)]\big{]} we obtain the correlation-picture generator as

[TABLE]

Although the dynamics described by $\mathpzc{D}_{\mathpzc{c}}$ is fully equivalent to the Schrödinger picture dynamics governed by $\mathpzc{D}_{\mathpzc{s}}$, working in the correlation picture offers clear advantages.

IV Universality of the Lindblad-like form for open system dynamics

We show in the following that working in the correlation picture leads to a universal Lindblad-like (ULL) equation. In the next sections, we shall discuss applications and further properties of the ULL formalism.

IV.1 General theory

From Eq. (6) we can readily obtain the dynamics of the subsystem by tracing over the bath degrees of freedom in $\dot{\varrho}_{\mathsf{SB}}(\tau)=\mathpzc{D}_{\mathpzc{c}}[\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)]$. To show that the subsystem dynamics has a Lindblad-like form, we use the expansions of $H_{\mathrm{I}}=\sum_{i=1}^{d_{\mathsf{S}}^{2}-1}\mathpzc{S}_{i}\otimes\mathpzc{B}_{i}$ and $H_{\chi}(\tau)=\sum_{j=0}^{d_{\mathsf{S}}^{2}-1}\mathpzc{S}_{j}\otimes\mathpzc{B}_{j}^{\chi}(\tau)$, where $\{\mathpzc{S}_{i}\}_{i=1}^{d_{\mathsf{S}}^{2}-1}$ is the basis of Hermitian operators defined on the system Hilbert space. Here, $d_{\mathsf{S}}$ is the dimension of the Hilbert space of the system and $\mathpzc{S}_{0}=\mathbbmss{I}$. We emphasize the difference between the operators $\mathpzc{B}_{i}$ and $\mathpzc{B}_{i}^{\chi}(\tau)$, which are related to $H_{\mathrm{I}}$ and $H_{\chi}$, respectively. Inserting these expansions into \dot{\varrho}_{\mathsf{S}}=\mathrm{Tr}_{\mathsf{B}}\big{[}-i[H_{\mathsf{SB}},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]-\big{[}H_{\mathsf{SB}},[\hskip-2.56073pt[H_{\chi},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]\hskip-2.56073pt]\big{]}\big{]} yields

[TABLE]

where

[TABLE]

are the elements of the covariance matrix $C$ of the bath operators (with $\langle{\circ}\rangle_{\mathsf{B}}=\mathrm{Tr}_{\mathsf{B}}[\varrho_{\mathsf{B}}(\tau)\circ]$). Here, unlike the standard Lindblad equation, these bath operators are obtained not only from the interaction Hamiltonian but also from the correlation parent operator $H_{\chi}$ defined in Eq. (3). We rewrite $C(\tau)$ in Eq. (7) in terms of its Hermitian and anti-Hermitian parts as $C(\tau)=A(\tau)+iB(\tau)$, with Hermitian matrices $A(\tau)=[a_{ij}(\tau)]$ defined by $a_{ij}(\tau)=(1/2)[C(\tau)+C^{{\dagger}}(\tau)]_{ij}$ and $B(\tau)=[b_{ij}(\tau)]$ defined by $b_{ij}(\tau)=(-i/2)[C(\tau)-C^{{\dagger}}(\tau)]_{ij}$, for $i,j\geqslant 1$. This leads to an exact Lindblad-like master equation for the system,

[TABLE]

where $\{X,Y\}\equiv XY+YX$ denotes the anticommutator. In this equation, the quasi-rates $\gamma_{m}^{\chi}$ are the eigenvalues of the matrix $A(\tau)$. The jump operators are given by $L_{m}^{\chi}=\sum_{j\neq 0}V_{mj}\mathpzc{S}_{j}$ where {$V_{mj}\}_{j}$ are the elements of the eigenvector corresponding to the eigenvalue $\gamma^{\chi}_{m}$. Furthermore, the Lamb-shift-like Hamiltonian {\mathbbmss{h}}_{{\mathsf{L}}}^{\chi}(\tau)=\langle{H_{\mathrm{I}}}\rangle_{\mathsf{B}}+2\textstyle{\sum_{i\neq 0}}\mathrm{Im}\big{(}c_{i0}(\tau)\big{)}\mathpzc{S}_{i}+\textstyle{\sum_{i\neq 0,j\neq 0}}b_{ij}(\tau)\mathpzc{S}_{i}\mathpzc{S}_{j}, where $b_{ij}$ are the elements of the matrix $B(\tau)$. For more details on the derivation see Ref. SM .

Several remarks are in order here:

(i) Although by combining known results in the literature one may infer time-local Lindblad-like forms for the dynamics under the assumptions of linearity or absence of initial system-environment correlations Andersson ; GKS ; book:Breuer-Petruccione , the ULL equation (9) is completely general; we have made no assumptions on the initial system-environment correlations or the strength of the system-environment interaction. In addition, the ULL equation is built explicitly on a microscopic theory of correlations in the total system.

(ii) Unlike the well-known Markovian embedding, where a non-Lindblad (non-Markovian) evolution is mapped to a Lindblad (Markovian) evolution for a specific larger system employing an ancillary system garrahan , we have proven here that the Lindblad-like form is derived for the open-system dynamics itself.

(iii) We note that the coefficients in the ULL master equation (8) refer to the correlation functions $\langle{\mathpzc{B}_{i}\mathpzc{B}_{j}^{\chi}(\tau)}\rangle_{\mathsf{B}}$, defined by the correlation parent operator and are thus different from the conventional correlation functions of bath operators $\langle{\mathpzc{B}_{i}(\tau)\mathpzc{B}_{j}(\tau)}\rangle_{\mathsf{B}}$ appearing in the standard Markovian Lindblad equation book:Breuer-Petruccione , where the average is taken on a constant bath state and $\mathpzc{B}_{i}(\tau)=U^{{\dagger}}_{\mathsf{B}}(\tau)\mathpzc{B}_{i}U_{\mathsf{B}}(\tau)$ [with $U_{\mathsf{B}}(\tau)=e^{-i\tau H_{\mathsf{B}}}$].

(iv) Since $\mathpzc{L}^{\chi}$ depends on the state of the system, Eq. (9) is formally a nonlinear equation. Indeed, the linearity constraint on the full dynamics of quantum systems does not imply a similar restriction on the dynamics of a subsystem, and this nonlinearity is naturally expected for a general dynamical equation. Nevertheless, we show in Ref. SM that our ULL master Eq. (9) is linear in two important cases: (a) if there is no initial correlation, i.e., $\chi(0)=0$, where we show $\chi(\tau)$ can be explicitly expressed in terms of the system-bath product state, and (b) if the domain of $\mathpzc{L}^{\chi}$ is restricted to a set of states $\{\varrho_{\mathsf{S}}^{(i)}\}$ forming a convex decomposition of the state of the system, i.e., $\varrho_{\mathsf{S}}=\sum_{i}p_{i}\varrho_{\mathsf{S}}^{(i)}$, but here the initial total state may be correlated.

IV.2 Example I: Jaynes–Cummings model with initial correlation

To illustrate universality of the dynamical Eq. (9), even in the presence of initial system–bath correlations, we begin with a proof-of-principle example, the well-known, exactly solvable Jaynes–Cummings model Book: Gerry-Knight , and show that the dynamics of the two-level system is described by the ULL equation even when the system is correlated with a bosonic mode.

Consider a two-level system interacting with a single bosonic mode under the Jaynes–Cummings Hamiltonian. The system Hamiltonian is $H_{\mathsf{S}}=(\omega_{0}/2)\sigma_{z}$, where $\sigma_{\pm}=(\sigma_{x}\pm i\sigma_{y})/2$ and $\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$ are the $x$, $y$, and $z$ Pauli operators, respectively, the bosonic bath Hamiltonian is $H_{\mathsf{B}}=\omega\hat{a}^{\dagger}\hat{a}$, where $\hat{a}^{\dagger}$ ($\hat{a}$) is the creation (annihilation) operator of the bosonic mode, and $H_{\mathrm{I}}=\lambda(\sigma_{+}\otimes\hat{a}+\sigma_{-}\otimes\hat{a}^{{\dagger}})$ describes the system–bath interaction. For simplicity, we assume that $\omega_{0}=\omega$ and that the initial state of the total system is in a correlated state $|\psi(0)\rangle=r_{1}|\mathrm{e},0\rangle+r_{2}|\mathrm{g},1\rangle$, where both $r_{1}$ and $r_{2}$ are real numbers. Choosing $\mathpzc{S}_{0}=\mathbbmss{I}/\sqrt{2}$, $\mathpzc{S}_{1}=\sigma_{x}/\sqrt{2}$, $\mathpzc{S}_{2}=\sigma_{y}/\sqrt{2}$, and $\mathpzc{S}_{3}=\sigma_{z}/\sqrt{2}$ as the basis operators, we can find $\{\mathpzc{B}^{\chi}_{i}\}_{i}$ (see Ref. SM ). Thus, the bath covariances are obtained as $c_{10}=c_{20}=0$, $c_{11}=c_{22}=\lambda(-2ir_{1}r_{2}+\alpha_{1}\alpha_{2})/(2\alpha_{1}^{2}-2)$, and $c_{12}=-c_{21}=\lambda(2r_{1}r_{2}\alpha_{1}+i\alpha_{2})/\left(2(1-\alpha_{1}^{2})\right)$, where $\alpha_{1}=(1-2r_{1}^{2})\cos(2\lambda\tau)$ and $\alpha_{2}=(1-2r_{1}^{2})\sin(2\lambda\tau)$. Following the steps of the derivation of Eq. (9), we obtain

[TABLE]

where \widetilde{\omega}_{0}=4\lambda r_{1}r_{2}\alpha_{1}/\big{(}1+4r_{1}^{2}-4r_{1}^{4}-(\alpha_{1}^{2}-\alpha_{2}^{2})\big{)}, \gamma_{1}^{\chi}=-\lambda\alpha_{2}/\big{(}2(1-\alpha_{1})\big{)}, and \gamma_{2}^{\chi}=\lambda\alpha_{2}/\big{(}2(1+\alpha_{1})\big{)} SM . We emphasize that this equation is in the ULL form and is valid even with initial system–bath correlations.

V Reduction to a Markovian equation

Based on our general dynamical equation where system–bath correlations are fully incorporated, we can obtain simpler expressions for the case where the correlations are small. This approach is valid, e.g., in the vicinity of time instants at which the correlation vanishes or becomes negligible. In other words, we introduce a weak-correlation approximation. In such cases, we can simplify our ULL master equation into a Markovian Lindblad-like (MLL) master equation, in which jump rates are positive—as expected from Markovian dynamics Maniscalco . Below, we show that this equation correctly characterizes the universal quadratic short-time behavior of the system dynamics where the standard Lindblad master equation may fail Rivas-short-time ; DelCampo .

We assume that at $\tau_{0}$ the correlation vanishes. Without loss of generality we take $\tau_{0}=0$, thus $\chi(0)=0$. We allow the correlations to accumulate in the subsequent time steps due to the dynamics. To first order in the time argument $\tau$, we find that the correlation satisfies Eq. (3) with $H_{\chi}(\tau)=\tau\widetilde{H}_{\mathrm{I}}(\tau)$, where $\widetilde{H}_{\mathrm{I}}(\tau)=\textstyle{\sum_{i\neq 0}}\mathpzc{S}_{i}\otimes(\mathpzc{B}_{i}-\langle{\mathpzc{B}_{i}}\rangle_{\mathsf{B}})-\textstyle{\sum_{i\neq 0}}\langle{\mathpzc{S_{i}}}\rangle_{\mathsf{S}}\mathpzc{B_{i}}$ and $\langle{\circ}\rangle_{\mathsf{S}}=\mathrm{Tr}_{\mathsf{S}}[\varrho_{\mathsf{S}}(\tau)\circ]$. Thus, from the knowledge of $H_{\chi}$, we can read $\mathpzc{B}_{j}^{\chi}(\tau)=\tau(\mathpzc{B}_{j}-\langle\mathpzc{B}_{j}\rangle_{\mathsf{B}})$, where $j\geqslant 1$. Substituting these expressions into Eq. (8) the bath covariance matrix becomes

[TABLE]

where $\mathrm{Cov}_{\mathsf{B}}(O_{1},O_{2})=\langle{O_{1}O_{2}}\rangle_{\mathsf{B}}-\langle{O_{1}}\rangle_{\mathsf{B}}\langle{O_{2}}\rangle_{\mathsf{B}}$. Since the covariance matrix $\bm{c}(\tau)$ is positive-semidefinite, $a_{ij}(\tau)=c_{ij}(\tau)$ and $b_{ij}(\tau)=0$. The positivity of $A$ implies positivity of the rates $\gamma_{m}^{\chi}\geqslant 0$, which is a necessary feature of a Markovian dynamical evolution. To obtain an equation with no dependence on the state of the bath (recall that $\mathbbmss{h}_{{\mathsf{L}}}^{\chi}(\tau)$ and $c_{ij}(\tau)$ depend on $\varrho_{\mathsf{B}}(\tau)$), we also expand $\varrho_{\mathsf{B}}(\tau)$ around $\tau_{0}=0$ and keep relevant terms up to the first order in $\tau$. Thus, we obtain

[TABLE]

where subscripts $\mathsf{B}_{0}$ and $\mathsf{S}_{0}$ indicate that the averages or covariances are taken with respect to $\varrho_{\mathsf{B}}(0)$ and $\varrho_{\mathsf{S}}(0)$ rather than $\varrho_{\mathsf{B}}(\tau)$ and $\varrho_{\mathsf{S}}(\tau)$. In Eq. (12), we have defined $\widetilde{H}_{\mathsf{B}}=H_{\mathsf{B}}+\langle{H_{\mathrm{I}}}\rangle_{\mathsf{S}_{0}}$ (see Ref. SM for more details). Equation (9)—bearing in mind Eq. (12)—describes the short-time Markovian dynamics around a point of vanishing correlation.

We emphasize that our weak-correlation assumption is exact up to the first order in $\tau$. If we extend this Markovian dynamical equation to longer times, it may still work as an approximation for the exact dynamics, e.g., when the correlation becomes repeatedly zero Lidar-Whaley ; SM . Although at first sight expanding around a point of vanishing correlation may seem equivalent to the standard Born approximation, we will illustrate in the next example that the MLL equation can be different from the Redfield equation. In addition, unlike the Redfield equation book:Breuer-Petruccione ; book:Rivas-Huelga , the MLL equation always keeps the state positive, hence avoiding the so-called slippage issue which afflicts the Redfield equation slippage ; Whitney .

A final remark regarding the applicability of the MLL approximation in other regimes is in order. If the system has an asymptotic state, a candidate for such a state can be $\varrho^{\star}=\sum_{n}\langle E_{n}|\varrho_{\mathsf{SB}}(0)|E_{n}\rangle|E_{n}\rangle\langle E_{n}|$, where $\{|E_{n}\rangle\}$ are the eigenvectors of $H_{\mathsf{SB}}$ SM . Now assume that (i) the system is strongly interacting with the bath, i.e., $H_{\mathrm{I}}$ is the dominant term in the total Hamiltonian such that $H_{\mathsf{SB}}\approx H_{\mathrm{I}}$, (ii) the interaction Hamiltonian contains only one term, $H_{\mathrm{I}}=\mathpzc{S}\otimes\mathpzc{B}$, and (iii) the initial state of the system and the bath is uncorrelated, $\chi(0)=0$. It is straightforward to show that under these conditions $\varrho^{\star}$ is an uncorrelated state, and thus, also in this case an MLL equation describes the asymptotic dynamics SM . In the next section we provide another case where the MLL approximation also holds, and later we illustrate these behaviors with two examples.

VI Dynamics of the correlation: systematic weak-correlation expansion

The above MLL approximation shows that the correlation picture and the ULL equation can offer far-reaching practical implications beyond their fundamental appeal. We identify that the weak-correlation approximation is the basic ingredient of the MLL equation. Expanding upon this is desirable as it can make the ULL methodology more amenable to practical investigations of a diverse set of systems.

VI.1 General theory

To go further and demonstrate that the ULL equation systematically enables such a rich approximative structure, in the following we develop a perturbative weak-correlation expansion for the ULL equation. In particular, by using the MLL toolbox as a starting point for a perturbative expansion of the correlation matrix in the interaction picture, we find an exact dynamical equation for the correlation $\bm{\chi}$ (boldface denotes the interaction picture) and expand it in terms of the interaction Hamiltonian $\bm{H}_{\mathrm{I}}$ as

[TABLE]

where $\bm{f}_{0}(\tau)\equiv\bm{\chi}(0)$ and $\bm{f}_{\ell}(\tau)$ is of the order of $\bm{H}^{\ell}_{\mathrm{I}}$, for $\ell\geqslant 1$. We incorporate the first $\ell+1$ terms within this expansion in the ULL dynamical equations for the system and the bath and derive the associated $\ell$th-order approximate ULL equations, referred to as the “ULL$\ell$” equations.

Let us provide an outline of the weak-correlation expansion—for details see Ref. SM . We start by the decomposition of the state of the total system as in Eq. (1). The total system at a later time $\tau+\xi$ is given by

[TABLE]

where $U_{\xi}=e^{-iH_{\mathsf{SB}}\xi}$. The first term of the above equation represents the evolution of an uncorrelated state $\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)$ which up to the first order in $\xi$ is given by the MLL dynamical equation. Using Eq. (129) and the definition of the correlation opertaor (1) at time $\tau+\xi$ one obtains after some algebra SM

[TABLE]

where we have defined

[TABLE]

Using integration and iterations a solution to Eq. (94) is obtained as

[TABLE]

which is symbolically in the form of Eq. (13). Starting from Eq. (96), we can systematically approximate the correlation and hence the ULL dynamical equation. If the initial system-bath correlation vanishes, the first-order approximation (with respect to $\widetilde{\bm{H}}_{\mathrm{I}}$) gives

[TABLE]

Using this relation to derive the correlation parent operator SM one can obtain a weak-correlation approximation for the ULL dynamical equation, which is of second order with respect to $\widetilde{\bm{H}}_{\mathrm{I}}$—hence it is referred to as the ULL$2$ equation. We need to solve the coupled differential equations for $\dot{\varrho}_{\mathsf{S}}(\tau)$ and $\dot{\varrho}_{\mathsf{B}}(\tau)$ in a self-consistent fashion, to obtain approximate states of the system and the bath—see Ref. SM for details and further elaborations.

For sufficiently short times Eq. (103) reduces to

[TABLE]

which is identical to the correlation operator obtained in the derivation of the MLL dynamical equation SM . Hence although in both ULL$2$ and the MLL equations $\chi$ is of the first order with respect to $\widetilde{\bm{H}}_{\mathrm{I}}$, the ULL$2$ equation can be expected to lead to an improvement over the performance of the MLL approximation SM .

Let us discuss another asymptotic-time case where the MLL approximation holds. This underlines that the utility of the MLL approximation is not necessarily limited to the short-time dynamics. If (i) the interaction is sufficiently weak (weak-coupling or weak-correlation regime), and (ii) the subsystem dynamics reaches a steady state in a finite time, one can approximate $\bm{\varrho}_{\mathsf{S}}(s)\otimes\bm{\varrho}_{\mathsf{B}}(s)$ in Eq. (103) with the tensor product of subsystem steady states $\bm{\varrho}_{\mathsf{S}}^{\star}\otimes\bm{\varrho}_{\mathsf{B}}^{\star}$ and hence $\widetilde{\bm{H}}_{\mathrm{I}}(s)\approx\widetilde{\bm{H}}_{\mathrm{I}}^{\star}$ at most times, which yields the MLL approximation $\bm{\chi}^{(1),\star}_{\mathrm{MLL}}(\tau)\approx-i\tau[\widetilde{\bm{H}}_{\mathrm{I}}^{\star},\bm{\varrho}_{\mathsf{S}}^{\star}\otimes\bm{\varrho}_{\mathsf{B}}^{\star}]$. Since Eq. (103) becomes exact in the weak-coupling limit and that typical quantum systems reach their steady state in finite times rapid-equilibration1 ; rapid-equilibration2 ; rapid-equilibration3 ; rapid-equilibration4 , one can conclude that the MLL equation modified by the subsystem steady states may hold at long times for typical systems SM .

Applying additional approximations on $\bm{\chi}^{(1)}$ by imposing time-locality and also assuming that the bath state remains constant in time, i.e., approximating $\bm{\chi}^{(1)}\approx-i\int_{0}^{\tau}ds\,[\widetilde{\bm{H}}_{\mathrm{I}}(s),\bm{\varrho}_{\mathsf{S}}(\tau)\otimes\bm{\varrho}_{\mathsf{B}}(0)]$, yields a time-local ULL$2$ equation.

Systematic derivation of higher-order correlation terms up to $\bm{\chi}^{(\ell-1)}$ yielding the ULL$\ell$ approximations is straightforward but algebraically heavy, and so we leave this for future work. For details see Ref. SM . However, similar to other approximate techniques, it usually suffices to consider the lowest order ULL$2$ or at most a few lowest orders. For a comparison with other techniques see also Ref. SM .

We emphasize that although the validity of the weak-correlation expansion hinges on the strength of the interaction Hamiltonian, it shows clear differences with standard weak-coupling approximations book:Breuer-Petruccione ; book:Rivas-Huelga . In particular, we note that our expansion uses the correlation in a direct fashion as the key ingredient.

In the following, we illustrate our weak-correlation ULL method through two examples. In example II we show that the MLL equation captures the exact dynamics with a good accuracy, and outperforms the standard Markovian Lindblad equation. In addition, the time-local ULL$2$ equation becomes tantamount to the second-order time-convolutionless (TCL$2$) dynamical equation for the system, which gives the exact dynamics for this example. We also compare our solutions with a coarse-graining (CG) method book:Rivas-Huelga ; Lidar-Whaley ; Rivas-short-time ; Schallerbook ; Schaller-Brandes ; Schaller-Brandes-2 ; Majenz , which fails to exceed the performance of the ULL$2$ solution. In example III, we show that $\bm{\chi}^{(1)}$ of Eq. (103) leads to a significant improvement in predicting the dynamics of the system whereas there the TCL$2$ and other approximating techniques such as the CG method do not provide such accuracy.

VI.2 Example II: Atom in a bosonic bath

Let us consider a two-level system (atom) interacting with a manymode bosonic bath initially in the thermal state $\varrho^{\beta}_{\mathsf{B}}=e^{-\beta\sum_{n}\omega_{n}\hat{a}^{{\dagger}}_{n}\hat{a}_{n}}/\mathrm{Tr}[e^{-\beta\sum_{n}\omega_{n}\hat{a}^{{\dagger}}_{n}\hat{a}_{n}}]$ at temperature $T=1/\beta$. Here $\hat{a}_{l}$ is the annihilation operator for mode $l$. The total Hamiltonian reads

[TABLE]

where $\mathpzc{O}_{\mathsf{B}}=\sum_{n}\kappa_{n}(\hat{a}_{n}+\hat{a}^{{\dagger}}_{n})$. Assuming that the atom at all times retains only a small correlation with the bath, we conclude that Eq. (12) applies and we obtain the following master equation:

[TABLE]

where $\gamma(\tau)=2\tau\,\mathrm{Cov}_{\mathsf{B}_{0}}(\mathpzc{O}_{\mathsf{B}},\mathpzc{O}_{\mathsf{B}})$, and \mathrm{Cov}_{\mathsf{B}_{0}}(\mathpzc{O}_{\mathsf{B}},\mathpzc{O}_{\mathsf{B}})=\textstyle{\int_{0}^{\infty}}d\omega\,J(\omega)\big{(}2n(\beta,\omega)+1\big{)} is given in terms of a spectral density function $J(\omega)$ and the bosonic occupation number $n(\beta,\omega)=(e^{\beta\omega}-1)^{-1}$. Equation (21) describes pure dephasing in the eigenbasis of $\sigma_{x}$ and gives the population of the excited state of the atom as

[TABLE]

The solution of the exact dynamics for this example has been provided in Ref. Braun for $\omega_{0}=0$ and under the assumption of an initial thermal state for the bath and an Ohmic spectral density for the couplings of the interaction Hamiltonian, $J(\omega)=\eta\omega(1+\omega^{2}/\omega^{2}_{c})^{-2}$, where $\omega_{c}$ is the cutoff frequency and $\eta$ denotes the coupling strength between the system and the bath. This provides a convenient means of studying the accuracy of Eq. (22). As we argued in Sec. V, under the mentioned conditions even in the highly strong-coupling regime, the MLL equation is exact in the asymptotic time. In this example, $H_{\mathrm{I}}$ has only one term, $H_{\mathsf{S}}$ is set to zero by $\omega_{0}=0$, and for the chosen spectral density and the uncorrelated initial state, all the conditions are satisfied. Hence the MLL equation gives an exact prediction for the asymptotic state, see Fig. 2.

Figure 2 shows the evolution of the excited-state population and compares our MLL and ULL$2$ methods with the Redfield equation (an equation obtained by applying only the weak-coupling and time-locality approximations on the exact dynamics book:Rivas-Huelga ), the TCL$2$ master equation book:Breuer-Petruccione ; book:Rivas-Huelga , and the exact solution. In addition, to make a comparison with a CG dynamical equation we used the results of Ref. Schaller-Brandes-2 ; see Fig. 2. For this particular example the time-local ULL$2$ dynamical equation is identical to the TCL$2$ dynamical equation SM , and both coincide with the exact dynamics. From this figure and the explicit form of the Redfield equation SM , it is clear that the Redfield equation is less accurate in the low-temperature limit. The MLL equation follows the exact solution relatively well, whereas the Redfield equation exhibits a relatively slower decay. Note that when $\omega_{0}=0$ the standard Lindblad equation is equivalent to the Redfield equation. For details of the derivation of the Redfield and the TCL$2$ equations and for the analysis of the short-time dynamics using the Lindblad-like model and the exact evolution, see Ref. SM .

In the following example, we illustrate that our Markovian approximation works well to describe the short-time dynamics in a system with non-Markovian features. We then show that the ULL$2$ equation outperforms other methods for later times and captures the long-time dynamics more accurately.

VI.3 Example III: Damped harmonic oscillator within a bath of oscillators

Consider a quantum harmonic oscillator interacting with a bath of oscillators with the total Hamiltonian given by

[TABLE]

where $M$ is the number of the bath oscillators—see Fig. 3. For simplicity of the analysis, we assume the initial system-bath state $(c_{0}|0\rangle+c_{1}|1\rangle)_{\mathsf{S}}\otimes|0\rangle^{\otimes M}_{\mathsf{B}}$, where $|i\rangle$ denotes the eigenstate of the corresponding number operator with eigenvalue $i$.

To obtain the MLL equation, we choose $\mathpzc{S}_{1}=\hat{a}+\hat{a}^{{\dagger}}$ and $\mathpzc{S}_{2}=i(\hat{a}-\hat{a}^{{\dagger}})$, and hence $\mathpzc{B}_{1}=\sum_{k}g_{k}(\hat{b}_{k}+\hat{b}_{k}^{{\dagger}})/2$ and $\mathpzc{B}_{2}=\sum_{k}ig_{k}(\hat{b}_{k}-\hat{b}_{k}^{{\dagger}})/2$. Inserting these into Eq. (12) yields $a_{11}=a_{22}=\tau G/4$, $a_{21}=a_{12}^{\ast}=i\tau G/4$, $\mathrm{Im}\langle{\mathpzc{B}_{1}\mathpzc{B}_{1}}\rangle_{\mathsf{B}_{0}}=\mathrm{Im}\langle{\mathpzc{B}_{2}\mathpzc{B}_{2}}\rangle_{\mathsf{B}_{0}}=0$, and $\mathrm{Im}\langle{\mathpzc{B}_{2}\mathpzc{B}_{1}}\rangle_{\mathsf{B}_{0}}=-\mathrm{Im}\langle{\mathpzc{B}_{1}\mathpzc{B}_{2}}\rangle_{\mathsf{B}_{0}}=G/4$, where $G=\sum_{k}|g_{k}|^{2}$; hence $\mathbbmss{h}_{\mathsf{L}}^{\chi}=0$. The MLL equation thus reads

[TABLE]

The standard Lindblad equation for this model has been given in Ref. Rivas-damped (see also Refs. Isar ; grabert for more general Lindblad forms for this model and their solutions). For the special case of the chosen initial state the Lindblad equation becomes

[TABLE]

where $\delta=\mathbbmss{P}[\int_{0}^{\infty}d\omega\,J(\omega)\,(\omega_{0}-\omega)^{-1}]$, $\gamma=\pi J(\omega_{0})$, $J(\omega)=\sum_{k=1}^{M}g_{k}^{2}\delta(\omega-\omega_{k})$ (the spectral density function), and $\mathbbmss{P}$ denotes the Cauchy principal value.

To compare the MLL equation with the standard Lindblad equation, the ULL$2$ equation, and the exact dynamics, we choose an Ohmic spectral density as $J(\omega)=(\omega/\pi)e^{-\omega/\omega_{c}}$ (similar to Ref. Rivas-damped ). For numerical simulation of this model, we take $\omega_{k}=0.1k$ and $M=255$. In order for the MLL equation and the simulations to be comparable with the standard Lindblad equation we choose the coupling strength as $g_{k}=\sqrt{J(\omega_{k})}$. It is shown in Fig. 4 (left) that the short-time dynamics is well captured by the MLL equation, for all chosen cutoff values, whereas the Lindblad equation fails to capture this and the TCL$2$ equation seems to capture only the initial moments of the dynamics. When $\omega_{c}=10$ the exact dynamics shows nondecaying oscillations, which may be due to non-Markovianity. It is interesting to note that even in this non-Markovian case the MLL equation—which is Markovian—can capture the short-time dynamics well. Considering the dynamics at long times, we observe in Fig. 4 (right) that the ULL$2$ equation can capture the exact dynamics with good accuracy, while the second-order NZ (NZ$2$) equation book:Breuer-Petruccione gives unphysical results. It is interesting that although the complete ULL and the NZ dynamical equations are exact, their second-order approximations, i.e., ULL$2$ and NZ$2$, approximate the dynamics differently. This is due to the different underlying approaches in applying the second-order approximation.

In addition, by calculating the population of the first excited state using the suggested asymptotic state $\varrho^{\star}$ SM , we obtain that $\varrho^{\star}_{11}=0.3963$ equals the asymptotic value in Fig. 4.

VII Summary and conclusions

We have introduced the correlation picture as a new dynamical picture. By using the Schrödinger equation of the total system, using a correlating transformation, and tracing over the environment degrees of freedom, we have found the dynamical equation of the subsystem without invoking any approximations. We have shown that this exact dynamical equation is in the Lindblad form, even if the system is initially correlated or is in the strong-coupling regime. Hence the Lindblad form for the dynamical equation is general, and the obtained master equation is a universal Lindblad-like (ULL) equation. We have provided a way to derive a Markovian master equation (MLL) from the ULL equation. In particular, we have shown that Markovianity can emerge if we apply a weak-correlation approximation, and the MLL equation becomes exact at instants with vanishing correlations. We have demonstrated that, not only at the initial time (if the system and the bath are prepared in a product state) but also in the asymptotic time, this weak-correlation approximation can be valid under certain conditions.

The correlation picture has also enabled us to formulate a systematic weak-correlation perturbative expansion, from which we have introduced approximate second- and higher-order tractable master equations. The MLL methodology plays an important role in this construction of the master equations which can feature non-Markovian effects. We have shown that existing and widely used weak-coupling-based equations emerge as special cases of our perturbative constructions, and thus our weak-correlation master equations are expected to outperform or perform as accurately as corresponding weak-coupling solutions. In particular, we have illustrated through three examples our results for the existence of the ULL equation, validity of the MLL equation around weak-correlation points for initial and asymptotic times, and have compared the MLL and ULL$2$ equations as approximate solutions with other Markovian and non-Markovian equations. We have shown that in these examples our equations describe the dynamics more accurately. We expect that introducing the correlation picture can pave the way for developing new techniques for controlling and harnessing system-environment correlations. We also anticipate a wide range of applications of our theory from quantum thermodynamics to quantum computation. In particular, our approach may help to understand whether and how quantum systems thermalize, and it may shed light on the role of correlations in quantum algorithms and the robustness of quantum error correction against correlated noise mechanisms.

Acknowledgements.

Discussions with J. Anders, E. Aurell, H.-P. Breuer, L. A. Correa, S. F. Huelga, A. Isar, S. J. Kazemi, J. Piilo, Á. Rivas, R. Sampaio, J. Tuorila, and S. Vinjanampathy are acknowledged. This work was supported by the Academy of Finland’s Center of Excellence program QTF Project 312298 and Sharif University of Technology’s Office of Vice President for Research and Technology. A.T.R. also acknowledges hospitality and support of the QTF Center of Excellence at Aalto University.

Appendix A Solution of the equation $\chi=-i[\hskip-2.56073pt[H_{\chi},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]\hskip-2.56073pt]$ for $H_{\chi}$

Equation (3) of the main text can be rewritten as

[TABLE]

which is of the general form $E^{{\dagger}}X+X^{{\dagger}}E=F$ (with unknown $X$), where we take $E=\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)$, $X=iH^{{\dagger}}_{\chi}(\tau)$, and $F=\chi(\tau)$. We define the instantaneous projection operator onto the null-space of $\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)$ by $P_{0}(\tau)$. From the result of Ref. [36] of the main text, if $P_{0}(\tau)\chi(\tau)P_{0}(\tau)=0$, then

[TABLE]

is a solution of Eq. (26), where $\varrho_{\mathsf{S}}^{-1}(\tau)$ and $\varrho_{\mathsf{B}}^{-1}(\tau)$ are the pseudo-inverses pseudo-invers of $\varrho_{\mathsf{S}}(\tau)$ and $\varrho_{\mathsf{B}}(\tau)$, respectively. The solution is not unique and may include further terms, \varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)Z(\tau)\Big{(}\mathbbmss{I}-P_{0}(\tau)\Big{)}+Y(\tau)P_{0}(\tau), where $Y(\tau)$ is an arbitrary operator and $Z(\tau)$ is an operator satisfying \varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)\big{(}Z(\tau)+Z^{{\dagger}}(\tau)\big{)}\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)=0. For our purposes in this paper, we take $Y(\tau)$ and $Z(\tau)$ equal to zero—with no impact on the state $\varrho_{\mathsf{SB}}$.

To show that Eq. (26) always has a solution, we need to ensure that the condition $P_{0}(\tau)\chi(\tau)P_{0}(\tau)=0$ is always satisfied, or equivalently, $P_{0}(\tau)\varrho_{\mathsf{SB}}(\tau)P_{0}(\tau)=0$, since by definition $P_{0}(\tau)\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)P_{0}(\tau)=0$. To this end we write $\varrho_{\mathsf{SB}}$ in terms of its spectral decomposition $\varrho_{\mathsf{SB}}=\sum_{i}\xi_{i}|\xi_{i}\rangle\langle\xi_{i}|$, and use the Schmidt decomposition (see Ref. [26] of the main text) for each eigenstate $\ket{\xi_{i}}=\textstyle{\sum_{j}}\sqrt{\lambda_{j}^{(i)}}\ket{e_{j}^{(i)}}_{\mathsf{S}}\ket{f_{j}^{(i)}}_{\mathsf{B}}$ to reach

[TABLE]

From $\mathrm{Tr}[P_{0}(\tau)\varrho_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)]=0$ one concludes that

[TABLE]

Since $\xi_{i}$ and $\lambda_{j}^{(i)}$ are positive numbers and $P_{0}$ is a positive–semidefinite matrix, for Eq. (29) to hold it is needed that $P_{0}(\tau)\ket{e_{j}^{(i)},f_{l}^{(k)}}_{\mathsf{SB}}=0$. Hence it is seen that applying $P_{0}(\tau)$ on Eq. (28) leads to $P_{0}(\tau)\varrho_{\mathsf{SB}}(\tau)P_{0}(\tau)=0$, which accordingly gives $P_{0}(\tau)\chi(\tau)P_{0}(\tau)=0$ .

Appendix B Derivation of the ULL equation

By tracing out over the bath from both sides of the dynamical equation in the correlation picture, we reach

[TABLE]

Inserting $H_{\mathsf{SB}}=H_{\mathsf{S}}+H_{\mathsf{B}}+H_{\mathrm{I}}$ and the identities (cyclicity of partial trace)

[TABLE]

into Eq. (30) yields

[TABLE]

where

[TABLE]

For future reference, we can also define $\widetilde{H}_{\mathsf{B}}(\tau)$ in a similar fashion. Expanding $H_{\mathrm{I}}$ in terms of a Hermitian operator basis $\{\mathpzc{S}_{i}\}_{i=1}^{d_{\mathsf{S}}^{2}-1}$ (with $\mathpzc{S}_{0}=\mathbbmss{I}$) as

[TABLE]

and replacing it into Eq. (32) yields

[TABLE]

We now expand $H_{\chi}(\tau)$ in the same basis

[TABLE]

and insert it into Eq. (35), which gives

[TABLE]

where it is seen that the elements of the matrix $C$ are given by $c_{ij}(\tau)=\mathrm{Tr}[\varrho_{\mathsf{B}}(\tau)\mathpzc{B}_{i}\mathpzc{B}_{j}^{\chi}(\tau)]$. Defining the Hermitian matrices $A=[a_{ij}]_{i\neq 0,j\neq 0}$ and $B=[b_{ij}]_{i\neq 0,j\neq 0}$, see the main text below Eq. (8), we obtain $c_{ij}=a_{ij}+i\,b_{ij}$ and $c^{\ast}_{ji}=a_{ij}-i\,b_{ij}$. Inserting these equations into Eq. (37) leads to the following general form for the dynamical equation:

[TABLE]

where

[TABLE]

The only remaining step to obtain the final form of Eq. (9) of the main text is to diagonalize $A=[a_{ij}]_{i\neq 0,j\neq 0}$, such that $VAV^{{\dagger}}=\Gamma$, where $\Gamma$ is a diagonal matrix with the eigenvalues $\gamma^{\chi}_{m}$ of $A$ as its diagonal elements.

Appendix C On the linearity of the ULL equation

Following the principles of quantum mechanics, the total dynamics should be linear for arbitrary initial states of the total system. However, linearity of the subsystem dynamics in general is not required. One can investigate the linearity of the Lindblad-like equation for the subsystem from different perspectives. Consider two reduced system states $\varrho^{(1)}_{\mathsf{S}}$ and $\varrho^{(2)}_{\mathsf{S}}$ obtained, respectively, by the partial trace over the bath of the total states $\varrho^{(1)}_{\mathsf{SB}}=\varrho^{(1)}_{\mathsf{S}}\otimes\varrho^{(1)}_{\mathsf{B}}+\chi^{(1)}$ and $\varrho^{(2)}_{\mathsf{SB}}=\varrho^{(2)}_{\mathsf{S}}\otimes\varrho^{(2)}_{\mathsf{B}}+\chi^{(2)}$. Now consider the convex combination of these states $p_{1}\varrho^{(1)}_{\mathsf{S}}+p_{2}\varrho^{(2)}_{\mathsf{S}}$, with $p_{1},p_{2}\geqslant 0$ and $p_{1}+p_{2}=1$, which is evidently obtained by the partial trace over the total system state $p_{1}\varrho^{(1)}_{\mathsf{SB}}+p_{2}\varrho^{(2)}_{\mathsf{SB}}$. Starting from the Schrödinger equation and following the derivation of the subsystem dynamics in Sec. B, it can be seen that in this case the generator of the Lindblad-like equation $\mathpzc{L}^{\chi}$ in general is not linear in the sense that

[TABLE]

where $\chi$ denotes the correlation operator when the total system state is $p_{1}\varrho^{(1)}_{\mathsf{SB}}+p_{2}\varrho^{(2)}_{\mathsf{SB}}$. However, below we discuss a restricted case where linearity in this sense can be retrieved.

C.1 Linearity on a restricted set of states

Consider a total state, defined with a given $\varrho_{\mathsf{B}}$ and $\chi$ such that $\varrho_{\mathsf{SB}}=\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}+\chi$, and initial subsystem states which are chosen from a restricted set of states forming a convex decomposition $\varrho_{\mathsf{S}}=\sum_{i}p_{i}\varrho_{\mathsf{S}}^{(i)}$, with $p_{i}\geqslant 0$ and $\sum_{i}p_{i}=1$. By replacing the convex decomposition of $\varrho_{\mathsf{S}}$ in Eq. (3) of the main text and defining $\chi^{(i)}:=-i[\hskip-2.56073pt[H_{\chi},\varrho_{\mathsf{S}}^{(i)}\otimes\varrho_{\mathsf{B}}]\hskip-2.56073pt]$, we observe that $\chi$ can also be written in the convex combination form of $\chi=\sum_{i}p_{i}\chi^{(i)}$. Thus, one can associate a correlation operator $\chi^{(i)}$ to each $\varrho^{(i)}_{\mathsf{S}}$ while $H_{\chi}$ remains identical for each of them. As a result, a single $\mathpzc{L}^{\chi}$ is associated with all of the cases. Following the derivation of the Lindblad-like equation and replacing convex decompositions of $\varrho_{\mathsf{S}}$ and $\chi$, one can obtain from the dynamical equation for $\varrho_{\mathsf{S}}$ the identity

[TABLE]

C.2 Linearity when there is no initial system-bath correlation

Now we consider another case where linearity holds. If the initial system-bath state is a product state, i.e., $\chi(0)=0$, and the dynamics of the total system is given by the unitary evolution $U(\tau)=e^{-iH_{\mathsf{SB}}}$, we obtain

[TABLE]

Thus, it follows that

[TABLE]

Replacing the above $\chi$ with $-i[\hskip-2.56073pt[H_{\chi},\varrho_{\mathsf{S}}\otimes\varrho_{\mathsf{B}}]\hskip-2.56073pt]$ in Eq. (30) yields

[TABLE]

It follows from the above equation that replacing $\varrho_{\mathsf{S}}(\tau)$ with $\alpha_{1}\varrho_{\mathsf{S}1}(\tau)+\alpha_{2}\varrho_{\mathsf{S}2}(\tau)$ and keeping $\varrho_{\mathsf{B}}(0)$ unchanged yield

[TABLE]

which implies that the reduced dynamics is linear.

Appendix D Details of example I: the Jaynes-Cummings model

Choosing $\mathpzc{S}_{0}=\mathbbmss{I}/\sqrt{2}$, $\mathpzc{S}_{1}=\sigma_{x}/\sqrt{2}$, $\mathpzc{S}_{2}=\sigma_{y}/\sqrt{2}$, $\mathpzc{S}_{3}=\sigma_{z}/\sqrt{2}$ as the system operator basis, we find $\mathpzc{B}_{0}=0$, $\mathpzc{B}_{1}=(\lambda/\sqrt{2})(\hat{a}+\hat{a}^{{\dagger}})$, $\mathpzc{B}_{2}=(i\lambda/\sqrt{2})(\hat{a}-\hat{a}^{{\dagger}})$, and $\mathpzc{B}_{3}=0$. Using the exact solution of the Jaynes-Cummings model (see Ref. [39] of the main text) we find |\psi(\tau)\rangle=e^{-i\tau\omega_{0}/2}\big{\{}\big{[}r_{1}\cos(\lambda\tau)-ir_{2}\sin(\lambda\tau)\big{]}|\mathrm{e},0\rangle+\big{[}-ir_{1}\sin(\lambda\tau)+r_{2}\cos(\lambda\tau)|\mathrm{g},1\rangle\big{]}\big{\}}, from which

[TABLE]

where $\alpha_{1}=(1-2r_{1}^{2})\cos(2\lambda\tau)$ and $\alpha_{2}=(1-2r_{1}^{2})\sin(2\lambda\tau)$. Now $H_{\chi}=\sum_{i=0}^{4}\mathpzc{S}_{i}\otimes\mathpzc{B}^{\chi}_{i}$, with

[TABLE]

Thus the bath covariances are obtained as $c_{10}=c_{20}=0$, $c_{11}=c_{22}=\lambda(-2ir_{1}r_{2}+\alpha_{1}\alpha_{2})/(2\alpha_{1}^{2}-2)$, and

[TABLE]

After obtaining $A$ and $B$ and diagonalizing $A$, we get

[TABLE]

Replacing these into Eq. (9) of the main text, the dynamical equation of the system is obtained in the ULL form.

Appendix E Derivation of the MLL equation

Using the definition of $\chi(\tau)$ from Eq. (1) of the main text and assuming $\chi(\tau_{0})=0$, we obtain

[TABLE]

from which, since $[\hskip-2.56073pt[H_{\chi}(\tau_{0}),\varrho_{\mathsf{S}}(\tau_{0})\otimes\varrho_{\mathsf{B}}(\tau_{0})]\hskip-2.56073pt]=i\chi(\tau_{0})=0$, we have

[TABLE]

and similarly for the bath,

[TABLE]

Using the above equations in $\dot{\chi}(\tau)=\dot{\varrho}_{\mathsf{SB}}(\tau)-\dot{\varrho}_{\mathsf{S}}(\tau)\otimes\varrho_{\mathsf{B}}(\tau)-\varrho_{\mathsf{S}}(\tau)\otimes\dot{\varrho}_{\mathsf{B}}(\tau)$ and expanding $\chi(\tau)$ around $\tau_{0}$ as $\chi(\tau_{0}+\tau)=\chi(\tau_{0})+\tau\dot{\chi}(\tau_{0})+O(\tau^{2})$ yields

[TABLE]

where (for an arbitrary $s$)

[TABLE]

where we remind that $\langle{\circ}\rangle_{\mathsf{S}}=\mathrm{Tr}[\varrho_{\mathsf{S}}(s)\circ]$ (and similarly for $\langle{\circ}\rangle_{\mathsf{B}}$). Equation (57) can also be written as

[TABLE]

which is more convenient for our analysis. Comparing this equation with Eq. (3) of the main text we conclude that

[TABLE]

Thus from Eq. (59) we conclude that for $j\neq 0$, $\mathpzc{B}_{j}^{\chi}=\tau(\mathpzc{B}_{j}-\langle\mathpzc{B}_{j}\rangle_{\mathsf{B}})$. Hence from Eq. (8) of the main text we obtain for $(i,j)\neq(0,0)$ that

[TABLE]

which is a positive matrix. Hence $A$ is positive and $B=0$. For $j=0$ we obtain $\mathpzc{B^{\chi}_{0}}=-\tau\textstyle{\sum_{i\neq 0}}\sqrt{d_{\mathsf{S}}}\langle{\mathpzc{S}_{i}}\rangle_{\mathsf{S}}\mathpzc{B}_{i}$, which yields

[TABLE]

and the dynamical equation is obtained as

[TABLE]

The above equation depends on the instantaneous state of the bath, which makes it not directly applicable as a system dynamical equation. To write it as an equation which depends only on the state of the system, we note that Eq. (64) is valid up to the second order in time around zero-correlation points. Thus, we can expand $\varrho_{\mathsf{B}}(\tau)$ around $\tau_{0}$ using Eq. (56), and keep only relevant terms. Replacing $\varrho_{\mathsf{B}}(\tau_{0}+\tau)=\varrho_{\mathsf{B}}(\tau_{0})-i\tau[\widetilde{H}_{\mathsf{B}}(\tau_{0}),\varrho_{\mathsf{B}}(\tau_{0})]+O(\tau^{2})$ and $\varrho_{\mathsf{S}}(\tau_{0}+\tau)=\varrho_{\mathsf{S}}(\tau_{0})-i\tau[\widetilde{H}_{\mathsf{S}}(\tau_{0}),\varrho_{\mathsf{S}}(\tau_{0})]+O(\tau^{2})$ into Eq. (62) yields

[TABLE]

where subscripts $\mathsf{S}_{0}$ and $\mathsf{B}_{0}$ mean that the averages are taken with respect to the states of the system and the bath at $\tau=\tau_{0}$. Replacing these into Eq. (37) leads to the following Markovian dynamical equation:

[TABLE]

where we have used

[TABLE]

Without loss of generality we can assume $\tau_{0}=0$.

Equation (66) has been obtained for sufficiently small $\tau$’s. However, if we extend this equation to arbitrary time instants, we shall have the MLL master equation, which is an approximation for the exact dynamics. This extension is a sort of coarse graining in time, yet with clear differences with the standard coarse-graining (CG) methods (see, e.g., Refs. [4, 17, 24, 46–49] of the main text). For example, the MLL coarse-graining and the CG method of Ref. [17] both assume that the correlation can vanish repeatedly during the evolution (a Born approximation), while MLL jump rates are proportional to the evolution time $\tau$—rather than a fixed CG time scale. In addition, a closer inspection shows that the way the Born approximation is implemented in our MLL equation is different from that of the CG methods (and also the Redfield equation). Methodologically, CG techniques often involve short-time, weak-coupling, and the standard Born approximation; whereas our MLL equation hinges on a weak-correlation approximation. Despite evident similarities, the final results appear different. We illustrate this point in the examples.

Appendix F The MLL equation is exact for short times

The short-time behavior of the system density matrix around $\tau_{0}$ is obtained by integration of Eq. (66), which yields

[TABLE]

To calculate $\textstyle{\int_{0}^{\tau}}\varrho(s)\,ds$ and $\textstyle{\int_{0}^{\tau}}s\varrho(s)\,ds$ for short times, we insert Eq. (68) into the integrals, and thus we obtain

[TABLE]

and

[TABLE]

Inserting Eqs. (69) and (70) into Eq. (68), the short-time dynamics of the system is obtained up to third order in time as

[TABLE]

where $\widetilde{H}_{\mathsf{S}}=H_{\mathsf{S}}+\langle{H_{\mathrm{I}}}\rangle_{\mathsf{B}_{0}}$ and $\mathbbmss{D}[\circ]=\sum_{i\neq 0,j\neq 0}\mathrm{Cov}(\mathpzc{B}_{i},\mathpzc{B}_{j})_{\mathsf{B}_{0}}(2\mathpzc{S}_{j}\circ\mathpzc{S}_{i}-\{\mathpzc{S}_{i}\mathpzc{S}_{j},\circ\})$. When the system-bath initial state is prepared in a product state (which is often the case), i.e., when $\chi(0)=0$, the short-time behavior of the dynamics, either Markovian or non-Markovian, will be given with the above equation. From this equation it is immediate that if the state of the system at $\tau_{0}$ commutes with $\widetilde{H}_{\mathsf{S}}(0)$, the system dynamics will be proportional to $\tau^{2}$. Otherwise the linear term in $\tau$ can dominate in short-time evolution. In the case where $[\widetilde{H}_{\mathsf{S}}(0),\varrho_{\mathsf{S}}(0)]=0$ and $[\widetilde{H}_{\mathsf{B}}(0),\varrho_{\mathsf{B}}(0)]=0$, Eq. (71) is simplified as

[TABLE]

Using the definition of $\mathrm{Cov}(\mathpzc{B}_{i},\mathpzc{B}_{j})_{\mathsf{B}_{0}}$ as given below Eq. (11) of the main text and the expansion of $H_{\mathrm{I}}$ from Eq. (34), after some straightforward algebra, the above equation can be recast as

[TABLE]

Appendix G Applicability of the MLL approximation in asymptotic regimes

Consider a quantum system which is coupled to an environment and the spectral decomposition of the total Hamiltonian is given by $H_{\mathsf{SB}}=\sum_{m}E_{m}|E_{m}\rangle\langle E_{m}|$. Let us define the time-averaged state of the total system as

[TABLE]

where we have assumed that the total Hamiltonian is time-independent. Since $\lim_{\tau\rightarrow\infty}(1/\tau)\textstyle{\int_{0}^{\tau}}ds~{}e^{i(\omega-\omega^{\prime})s}\approx\delta_{\omega,\omega^{\prime}}$ and with the assumption that there is no degeneracy in the total Hamiltonian, it is seen that the asymptotic limit of $\bar{\varrho}_{\mathsf{SB}}(\tau)$ is obtained as

[TABLE]

From this equation it is evident that $[H_{\mathsf{SB}},\varrho_{\mathsf{SB}}^{\star}]=0$, hence $\varrho_{\mathsf{SB}}^{\star}$ is a steady state of the total system. Noting that, for generic system-bath Hamiltonians (with nondegenerate energies and gaps), for almost all observables $A$, we have $\langle A\rangle_{\varrho_{\mathsf{SB}}(\infty)}\approx\langle A\rangle_{\varrho^{\star}_{\mathsf{SB}}}$ Reimann , in this sense one may consider the state $\varrho^{\star}_{\mathsf{SB}}$ effectively as a description for the asymptotic state $\varrho_{\mathsf{SB}}(\infty)$. Hence the asymptotic state of the subsystem $\mathsf{S}$ can be read as $\varrho_{\mathsf{S}}^{\star}=\sum_{m}\langle E_{m}|\varrho_{\mathsf{SB}}(0)|E_{m}\rangle\mathrm{Tr}_{\mathsf{B}}{\left[|E_{m}\rangle\langle E_{m}|\right]}$.

In the highly strong-coupling regime, the total Hamiltonian can be approximated as $H_{\mathsf{SB}}\approx H_{\mathrm{I}}$. If the interaction Hamiltonian has only one term, i.e., $H_{\mathrm{I}}=\mathpzc{S}\otimes\mathpzc{B}$, the eigenbasis of the total Hamiltonian can be approximated with the eigenbasis of the interaction Hamiltonian. By decomposing $\mathpzc{S}$ and $\mathpzc{B}$ as $\mathpzc{S}=\sum_{i}\alpha_{i}|\alpha_{i}\rangle\langle\alpha_{i}|$ and $\mathpzc{B}=\sum_{j}\eta_{j}|\eta_{j}\rangle\langle\eta_{j}|$, the eigenvector of the total Hamiltonian $|E_{m}\rangle$’s can be approximated as $|E_{ij}\rangle=|\alpha_{i}\rangle\otimes|\eta_{j}\rangle$. Now by assuming that the initial state is uncorrelated, it is readily seen that $\varrho_{\mathsf{SB}}^{\star}$ is uncorrelated,

[TABLE]

Hence the correlation operator in this long-time limit vanishes.

Remark 1

As is clear from the above analysis, the MLL equation for long times can be generally different from the MLL equation obtained by initial conditions. In addition, constructing the MLL equation for long times requires the knowledge of the asymptotic state. However, this information is initially unavailable and its determination is indeed a goal of any dynamical equation in the first place. Despite this difficulty, for some systems the very MLL equation obtained by initial conditions may still provide the asymptotic state $\varrho^{\star}_{\mathsf{S}}$ with good accuracy. Example II in Sec. I provides a case of this type [see also Fig. 2 (left) in the main text]. In Sec. H.3.1 we show another time-asymptotic regime in which the MLL approximation still holds, but in the weak-correlation limit.

Appendix H Details of the derivation of the exact dynamical equation for the correlation operator and the systematic weak-correlation expansion

H.1 Exact dynamical equation for the correlation

Using the decomposition of the state of the total system at a given time $\tau_{0}$ as $\varrho_{\mathsf{SB}}(\tau_{0})=\varrho_{\mathsf{S}}(\tau_{0})\otimes\varrho_{\mathsf{B}}(\tau_{0})+\chi(\tau_{0})$, the state of the total system at a later time $\tau_{0}+\xi$ is given by

[TABLE]

where $U_{\xi}$ is the evolution operator of the total system $\mathsf{SB}$ in a time interval $\xi$. In the following we use the notation $\hat{\circ}$ to indicate that the related operator has been obtained from an uncorrelated state—which is akin to the MLL method. Let us remind that in our MLL methodology (Sec. E) we have already shown that the correlation formed between the system and the bath due to the evolution from an uncorrelated state (here $\varrho_{\mathsf{S}}(\tau_{0})\otimes\varrho_{\mathsf{B}}(\tau_{0})$) is given by [Eq. (57)]

[TABLE]

which is exact for short times $\xi$’s. Since the exact dynamical equation for the system and the bath are given by

[TABLE]

if we start from a zero-correlation state at $\tau_{0}$, at a later time $\tau_{0}+\xi$ the system and bath states are given up to the first order in time $\xi$ by

[TABLE]

Thus

[TABLE]

where on the right-hand side (RHS) we should keep terms of the first order in $\xi$. From this relation the state of the system (bath) is obtained by tracing out over the bath (system) as

[TABLE]

From Eqs. (83) – (85) and using the relation $\chi(\tau+\tau_{0})=\varrho_{\mathsf{SB}}(\tau_{0}+\xi)-\varrho_{\mathsf{S}}(\tau_{0}+\xi)\otimes\varrho_{\mathsf{B}}(\tau_{0}+\xi)$ the system-bath correlation operator at $\tau_{0}+\xi$ (up to the first order in $\xi$) is obtained as follows:

[TABLE]

or equivalently,

[TABLE]

Equations (84), (85), and (87) indicate that starting at time $\tau_{0}$ with the knowledge of the total state $\varrho_{\mathsf{SB}}(\tau_{0})$ (hence $\varrho_{\mathsf{S}}(\tau_{0})$, $\varrho_{\mathsf{B}}(\tau_{0})$, and $\chi(\tau_{0})$) one can read the corresponding quantities at a close later time $\tau_{0}+\xi$. When $\chi(\tau_{0})=0$, the states $\varrho_{\mathsf{S}}(\tau_{0}+\xi)$ and $\varrho_{\mathsf{B}}(\tau_{0}+\xi)$ only require the knowledge of the system and the bath at time $\tau_{0}$. However, the correlation operator $\chi(\tau_{0}+\xi)$ is required for the estimation of the states $\varrho_{\mathsf{S}}(\tau_{0}+2\xi)$ and $\varrho_{\mathsf{B}}(\tau_{0}+2\xi)$. Figure 5 depicts this hierarchical and algorithmic construction.

One can rewrite the above equation as a recursive formula for calculating $\chi$ at a given time $\tau_{n+1}$ through an expression for time $\tau_{n}$ as

[TABLE]

From the above equation and by taking the $\xi\rightarrow 0$ limit, one can also derive the exact and general formula for $\dot{\chi}(\tau)$,

[TABLE]

To simplify the notation, we introduce a linear superoperator $\mathcal{Y}_{\tau}$ as

[TABLE]

using which Eq. (89) can be rewritten as

[TABLE]

This exact equation can be solved formally by iteration,

[TABLE]

where $s_{0}\equiv\tau$, $k$ counts the number of nested integrals, and the $k=0$ terms in the summations simply denote the identity map. This is an exact relation which shows how by having the exact subsystem states $\varrho_{\mathsf{S}}(\tau)$ and $\varrho_{\mathsf{B}}(\tau)$, the interaction Hamiltonian $H_{\mathrm{I}}$, and the initial correlation $\chi(0)$ one can read the development of the correlation operator with time.

We remark that an alternative way to derive Eq. (91) is to start from the very definition of the correlation operator $\chi(\tau)$ [Eq. (1) of the main text], differentiating it with respect to $\tau$ and using Eqs. (79) and (80), which yield Eq. (91). However, a clear physical appeal of the recursive-relation approach we developed above is that this way one can inspect how the development of correlation features beyond the Markovian approximation (MLL) play a key role in deriving more accurate dynamical equation. In addition, this approach better manifests the intricate role of the MLL expansion in obtaining an expansion for $\chi(\tau)$.

Two limiting cases can be discerned: (i) when there is no coupling between the system and the bath ($H_{\mathrm{I}}=\widetilde{H}_{\mathrm{I}}(\tau)=0$), this relation gives $\chi(\tau)=U_{\mathsf{S}}(\tau)\otimes U_{\mathsf{B}}(\tau)\chi(0)U_{\mathsf{S}}^{{\dagger}}(\tau)\otimes U^{{\dagger}}_{\mathsf{B}}(\tau)$, which obviously vanishes for the uncorrelated initial state. (ii) When there is no initial correlation $\chi(0)=0$ (but $H_{\mathrm{I}}\neq 0$ ), Eq. (92) reduces to the solution reported in Eq. (44).

As we delineate in the next section, the principal importance of Eq. (92) for our framework is that it also enables a systematic approach to approximating $\chi$ perturbatively—and hence the ULL master equation. To rigorously apply a weak-correlation expansion, it is important to derive the interaction-picture forms of the above equations.

H.2 Systematic weak-correlation perturbative expansion

The interaction picture forms of Eqs. (88) – (92) are as follows (we use boldfaced letters for operators in the interaction picture):

[TABLE]

where

[TABLE]

with $U_{0}(s)=e^{-is(H_{\mathsf{S}}+H_{\mathsf{B}})}=U_{\mathsf{S}}(s)\otimes U_{\mathsf{B}}(s)$. We note that in the first summation of Eq. (96) the $k$th term is of the order of $\bm{H}^{k}_{\mathrm{I}}$, whereas in the second summation the $k$th term is of the order of $\bm{H}^{k+1}_{\mathrm{I}}$. In addition, we can write solution (96) as

[TABLE]

where $\bm{f}_{0}(\tau)\equiv\bm{\chi}(0)$ and $\bm{f}_{\ell}(\tau)$ (for $\ell\geqslant 1$) is the addition of the $\ell$th term of the first summation and the $(\ell-1)$th term of the second summation in Eq. (96) [hence $\bm{f}_{\ell}$ is of the order of $\bm{H}^{\ell}_{\mathrm{I}}$]. This series solution to Eq. (94) is convergent as long as the kernel of its corresponding integral equation (the homogenous part) is square integrable book:Arfken . This can be readily satisfied for systems with $\|\bm{\mathcal{Y}}_{\tau}\|<\infty$—or simply $\|\bm{H}_{\mathrm{I}}\|<\infty$.

H.2.1 ULL$2$ equation

From Eq. (96) we can obtain the correlation operator up to the first order in the interaction Hamiltonian $\widetilde{\bm{H}}_{\mathrm{I}}$ as

[TABLE]

For an uncorrelated initial state ($\bm{\chi}(0)=0$) the above equation reduces to

[TABLE]

Inserting this relation in the interaction picture form of the dynamical equation (79) yields the following second-order equation:

[TABLE]

Returning to the Schrödinger picture, we have

[TABLE]

Inserting this into the system dynamical equation (79) gives

[TABLE]

which can also be written equivalently as

[TABLE]

We need to solve this equation together with the similar equation for $\dot{\varrho}_{\mathsf{B}}(\tau)$,

[TABLE]

We call the dynamical equations (106) [or Eq. (104)] and (108) the “ULL$2$” equations since they have terms of the order of $\bm{H}^{2}_{\mathrm{I}}$ (although they have been obtained from the first-order weak-correlation approximation $\bm{\chi}^{(1)}$).

Remark 2

In the case of a correlated initial state the analysis is similar to that above with extra terms in the relations.

Remark 3

From the expression for $\bm{\chi}^{(1)}$, one can also obtain an approximate expression for the correlation generator $\bm{H}^{(1)}_{\chi}$ by inserting Eq. (103) into the interaction-picture form of Eq. (27).

H.2.2 Higher orders

Higher order approximations to the ULL equation (referred to as “ULL$k$”) can be obtained by using $\bm{\chi}^{(k-1)}$ (for $k\geqslant 3$) given by Eq. (96). Deriving the explicit form of the ULL$k$ equation is straightforward, but the expressions become considerably involved. Fortunately, for many problems of interest it may suffice to consider the ULL$2$ or simply a few lowest orders, and we leave a more comprehensive analysis of the higher order ULL$k$ equations for future investigations.

H.3 Comparing the ULL$2$ equation and corresponding weak-coupling equations

There exist a number of powerful perturbative techniques in the literature for open systems. Here we discuss how some of these methods can be compared to or regarded as special cases of the ULL$2$ equation. In particular, we discuss the MLL equation, the second-order Nakajima-Zwanzig (NZ$2$) equation, and the second-order time-convolutionless (TCL$2$) equation. In addition, we also compare the ULL$2$ equation with the coarse-grained (CG) master equations Schaller-CG ; Majenz for the two examples discussed in the main text—and also here in Secs. I and J. A more comprehensive and detailed comparative analysis with other existing methods is beyond the scope of this paper.

H.3.1 Comparing the ULL$2$ equation with the MLL equation: ULL$2$ as an enhanced MLL

For short times Eq. (103) reduces to

[TABLE]

which is exactly the correlation operator obtained in the derivation of the MLL dynamical equation [cf. Eq. (60)]. Hence although in both ULL$2$ and the MLL equations $\chi$ is of the first order with respect to $\widetilde{\bm{H}}_{\mathrm{I}}$ [see Eqs. (103) and (57)], the ULL$2$ equation can expectedly lead to an enhancement over the performance of the MLL approximation and it can also include non-Markovian effects. In this sense, the ULL$2$ equation is an enhanced MLL dynamical equation.

As is evident from the construction of the MLL approximation, and also recalling Sec. G, the utility of the MLL approximation is not necessarily limited to the short-time dynamics. Here we discuss yet another case where the MLL approximation can be applied. If the following conditions hold: (i) the weak-correlation or weak-coupling assumptions hold (in contrast to the strong-coupling case of Sec. G), and (ii) the subsystem dynamics reaches a steady state in a finite time, one can approximate $\bm{\varrho}_{\mathsf{S}}(s)\otimes\bm{\varrho}_{\mathsf{B}}(s)$ in Eq. (103) with the tensor product of subsystem steady states $\bm{\varrho}_{\mathsf{S}}^{\star}\otimes\bm{\varrho}_{\mathsf{B}}^{\star}$ and hence $\widetilde{\bm{H}}_{\mathrm{I}}(s)\approx\widetilde{\bm{H}}_{\mathrm{I}}|_{\varrho_{\mathsf{S}}(s)\approx\varrho^{\star}_{\mathsf{S}},\varrho_{\mathsf{B}}(s)\approx\varrho^{\star}_{\mathsf{B}}}=:\widetilde{\bm{H}}^{\star}_{\mathrm{I}}$ at most times, which yields an MLL approximation $\bm{\chi}^{(1)\,\star}_{\mathrm{MLL}}(\tau)\approx-i\tau[\widetilde{\bm{H}}_{\mathrm{I}}^{\star},\bm{\varrho}_{\mathsf{S}}^{\star}\otimes\bm{\varrho}_{\mathsf{B}}^{\star}]$. Since Eq. (103) is accurate in the weak-correlation (or weak-coupling) regime, and typical quantum open systems often reach their steady state in finite times (see Refs. [42–45] of the main text), one can conclude that the MLL equation modified by the steady state still holds for long times for such typical systems.

Remark 4

Note that the same issue as in Remark 1 applies here too; a priori knowledge of the subsystem steady state seems necessary to construct the associated MLL equation. However, as shall be discussed later in Sec. H.4, one may alleviate this issue by applying the MLL approximation for both system and bath states around the initial time and solving these coupled equations simultaneously. Figure 6 shows that these initial-time MLL equations can give a significantly more accurate asymptotic state than other techniques (compare also with Fig. 4 of the main text where we report the result of other techniques).

H.3.2 Comparing the ULL$2$ equation with the NZ$2$ equation

We recall that the NZ$2$ equation (see, e.g., Ref. [3] of the main text) reads as

[TABLE]

under the assumption of $\mathrm{Tr}_{\mathsf{B}}[\bm{H}_{\mathrm{I}}(s)\bm{\varrho}_{\mathsf{B}}(0)]=0$. We observe that the ULL$2$ equation [Eq. (104)] is to some extent similar to the NZ$2$ equation. But two key differences can also be discerned by comparing the integrands of Eqs. (104) and (110): (i) the state of the bath is time-dependent in the ULL$2$ equation, whereas in the NZ$2$ equation it is a constant state, which is usually assumed to be the initial (or equilibrium) state of the bath, and (ii) in the ULL$2$ equation rather than $\bm{H}_{\mathrm{I}}(s)$ we have its effective form $\widetilde{\bm{H}}_{\mathrm{I}}(s)$. These points suggest that the performance of the ULL$2$ equation is generally different from that of the NZ$2$ equation, and that one may expect that the ULL$2$ equation either performs similarly to the NZ$2$ equation or outperforms it. To illustrate the validity of this prediction, we compare these two equations for example III of the main text (Sec. J of this supplementary material)—in particular, see Fig. 4 of the main text. There we show that the NZ$2$ equation may even give unphysical solutions whereas the ULL$2$ equation gives fairly accurate results.

H.3.3 Comparing the ULL$2$ equation with the TCL$2$ equations

We recall that the TCL$2$ equation for a general system in the interaction picture is (see, e.g., Refs. [3,4] of the main text)

[TABLE]

under the assumption of $\mathrm{Tr}_{\mathsf{B}}[\bm{H}_{\mathrm{I}}(s)\bm{\varrho}_{\mathsf{B}}(0)]=0$. Basically this equation has been obtained by a time-local transformation in the steps of the NZ framework—compare Eqs. (111) and (110). In order to make the ULL$2$ equation comparable (104) with the TCL$2$ equation, we first need to introduce a time-local version of the ULL$2$ equation.

H.3.4 Time-local ULL$2$ equation

Here we argue that by applying appropriate approximations on the ULL$2$ equation (104), one can obtain a dynamical equation which is comparable to the TCL$2$ equation. Let us start from the correlation operator $\bm{\chi}^{(1)}(\tau)$ [Eq. (103)]. If we apply the time-locality (or, in some sense, Markovian) approximation $\bm{\varrho}_{\mathsf{S}}(s)\otimes\bm{\varrho}_{\mathsf{B}}(s)\approx\bm{\varrho}_{\mathsf{S}}(\tau)\otimes\bm{\varrho}_{\mathsf{B}}(\tau)$ on the RHS of Eq. (103) and insert this correlation in the ULL$2$ equation, we obtain the following time-local ULL$2$ equation:

[TABLE]

If we further assume (as approximation) that the state of the bath does not change appreciably in time, ${\bm{\varrho}}_{\mathsf{B}}(\tau)\approx{\bm{\varrho}}_{\mathsf{B}}(0)$, and also that $\mathrm{Tr}_{\mathsf{B}}[\bm{H}_{\mathrm{I}}(\tau)\bm{\varrho}_{\mathsf{B}}(0)]=0$ (similar to the derivation of the TCL equation book:Rivas-Huelga ), the above time-local dynamical equation reduces to

[TABLE]

This equation is to some extent similar to the TCL$2$ equation (111)—but not fully. Due to the nature of the approximations made to derive this reduced time-local ULL$2$ equation, it is expected that the ULL$2$ equation in general may outperform the TCL$2$ equation, except under particular circumstances where the TCL$2$ equation may become exact.

Remark 5

Let us start from the correlation operator $\bm{\chi}^{(1)}(\tau)$ [Eq. (103)] and compare it with the definition of $\bm{\chi}(\tau)=-i[\bm{H}_{\chi}(\tau),\bm{\varrho}_{\mathsf{S}}(\tau)\otimes\bm{\varrho}_{\mathsf{B}}(\tau)]$ [Eq. (3) of the main text in the interaction picture]. If we apply the time-locality approximation by replacing $\bm{\varrho}_{\mathsf{S}}(s)\otimes\bm{\varrho}_{\mathsf{B}}(s)$ in the commutator of the integrand of Eq. (103) with its time-local form $\bm{\varrho}_{\mathsf{S}}(\tau)\otimes\bm{\varrho}_{\mathsf{B}}(\tau)$, we can read the correlation generator as $\bm{H}_{\chi}^{(1)\,\mathrm{TL}}(\tau)=\textstyle{\int_{0}^{\tau}}ds\,\widetilde{\bm{H}}_{\mathrm{I}}(s)$. If we also apply the time-locality approximation in $\widetilde{\bm{H}}_{\mathrm{I}}(s)$, the integrand in the above relation changes to $\widetilde{\bm{H}}_{\mathrm{I}}(s)|_{\bm{\varrho}_{\mathsf{S}}(s)\approx\bm{\varrho}_{\mathsf{S}}(\tau),\bm{\varrho}_{\mathsf{B}}(s)\approx\bm{\varrho}_{\mathsf{B}}(\tau)}$ [see Eq. (58) in the interaction picture]. Returning to the Schrödinger picture, we obtain $H_{\chi}^{(1)\,\mathrm{TL}}(\tau)=\textstyle{\int_{0}^{\tau}}ds\,U_{0}(\tau-s)\widetilde{H}_{\mathrm{I}}(s)|_{\varrho_{\mathsf{S}}(s)\approx\varrho_{\mathsf{S}}(\tau),\varrho_{\mathsf{B}}(s)\approx\varrho_{\mathsf{B}}(\tau)}U_{0}(s-\tau)$. Consequently, if we insert the above relation in the dynamical equation of the system [Eq. (32)] we obtain

[TABLE]

H.4 The dependence of the dynamical equation of the system on the instantaneous state of the bath

It is clear from Eq. (96) that the correlation operator depends on the instantaneous state of the bath. Hence the correlation generator $H_{\chi}$ and $\dot{\varrho}_{\mathsf{S}}(\tau)$ also explicitly depend on $\varrho_{\mathsf{B}}(\tau)$, and the ULL equation is more general but also more complicated than the NZ and TCL equations. Since our weak-correlation methodology hinges on obtaining or approximating the correlation operator (which is a joint property of the system and the bath), once we have an expression for $\chi$, we can employ it to also approximate the state of the bath more accurately. We can then use this approximate state of the bath to obtain a yet more accurate approximation for the dynamical equation of the system. Although this increases the number of the equations to be solved, it enables a more accurate estimation of the system. We will explicitly demonstrate this improvement in the next sections—Secs. I and J—where we elaborate on examples II and III of the main text.

We also show that it may not always be even approximately adequate to assume the bath to remain unchanged ($\varrho_{\mathsf{B}}(\tau)\approx\varrho_{\mathsf{B}}(0)$) during the dynamics of the system. In particular, in Fig. 6, we observe that this assumption leads to unphysical, negative populations in our example III, whereas using our MLL equation for the bath state too provides a considerable improvement over the constant-bath assumption. The situation can be improved further by using the ULL$2$ equations for both the system and the bath—Fig. 8 (and also Fig. 4 of the main text). We remark that other approximate methods may also provide some improvements over the constant-bath assumption.

Appendix I Details of example II: Two-level system in a bosonic bath

Here we provide the details related to example II introduced in Sec. VI A of the main text.

I.1 Redfield equation

Following the standard steps for derivation of the Redfield equation (e.g., see Ref. [3] of the main text) we obtain for the dynamics of a two-level system (atom)

[TABLE]

where $S(\beta,\omega)=2[n(\beta,\omega)+1]J(\omega)$. This yields the following dynamical equations for the populations:

[TABLE]

The solutions of these equations are given by

[TABLE]

To compare the Redfield solution with the exact solution, we rewrite the Redfield equation in the $\omega_{0}=0$ limit, which becomes equal to the Markovian Lindblad equation. Using the relation $\lim_{\omega_{0}\rightarrow 0}S(\beta,\omega_{0})=2\eta/\beta$, Eq. (115) reduces to

[TABLE]

Thus, we can obtain the differential equations for the diagonal elements of the density matrix,

[TABLE]

where $l\in\{\mathrm{g},\mathrm{e}\}$. The solutions of these Redfield equations are given by

[TABLE]

I.2 Time-local ULL$2$ equation

If in the time-local ULL$2$ equation (113) we calculate the averages in $\widetilde{\bm{H}}_{\mathrm{I}}(\tau)$ with respect to the initial system-bath state $\bm{\varrho}_{\mathsf{SB}}(0)=|\mathrm{e}\rangle\langle\mathrm{e}|\otimes\varrho^{\beta}_{\mathsf{B}}$ and put $\langle\circ\rangle_{\mathsf{S}(\mathsf{B})}\approx\langle\circ\rangle_{\mathsf{S}_{0}(\mathsf{B}_{0})}$, we arrive at $\widetilde{\bm{H}}_{\mathrm{I}}(\tau)\approx{\bm{H}}_{\mathrm{I}}(\tau)$, $\widetilde{\bm{H}}_{\mathsf{S}}(\tau)\approx{\bm{H}}_{\mathsf{S}}(\tau)$, which reduces the time-local ULL$2$ equation to

[TABLE]

This is identical to the TCL$2$ equation (111)—which is discussed in detail below.

I.3 TCL$2$ equation

Following the derivation of the TCL$2$ equation (e.g., see Ref. [4] of the main text) we obtain

[TABLE]

where boldface denotes the interaction picture, e.g., $\bm{\sigma}_{x}(s)=U^{{\dagger}}_{0}(s)\sigma_{x}U_{0}(s)$, $G(\tau,\tau-u)=\langle{\bm{\mathcal{O}}_{\mathsf{B}}(\tau)\bm{\mathcal{O}}_{\mathsf{B}}(\tau-u)}\rangle_{{\mathsf{B}}_{0}}$, and $G(\tau-u,\tau)=\langle{\bm{\mathcal{O}}_{\mathsf{B}}(\tau-u)\bm{\mathcal{O}}_{\mathsf{B}}(\tau)}\rangle_{{\mathsf{B}}_{0}}$. Since $\omega_{0}=0$, the interaction picture and Schrödinger picture are equivalent for the system operators; hence the above equation can be rewritten as

[TABLE]

Noting that the bath initial state is a thermal state and using a spectral density $J(\omega)$, we have G(\tau,\tau-u)=\int_{0}^{\infty}d\omega\,J(\omega)\,[\big{(}n(\beta,\omega)+1\big{)}e^{-i\omega\tau}+n(\beta,\omega)e^{i\omega\tau}] and G(\tau-u,\tau)=\int_{0}^{\infty}d\omega\,J(\omega)[\big{(}n(\beta,\omega)+1\big{)}e^{i\omega\tau}+n(\beta,\omega)e^{-i\omega\tau}], where $n(\beta,\omega)=1/(e^{\beta\omega}-1)$. Using these quantities the TCL$2$ dynamical equation is given by

[TABLE]

where $\gamma(\tau)=2\textstyle{\int_{0}^{\tau}}d\tau\,\int_{0}^{\infty}d\omega\,J(\omega)\cos(\omega\tau)[2n(\beta,\omega)+1]$.

I.4 Short-time dynamics

Since $\langle{H_{\mathrm{I}}}\rangle_{\mathsf{B}_{0}}=0$ we have $\widetilde{H}_{\mathsf{S}}=H_{\mathsf{S}}$. Thus, the initial state of the atom $\varrho_{\mathsf{S}}(0)=|\mathrm{e}\rangle\langle\mathrm{e}|$ commutes with $\widetilde{H}_{\mathsf{S}}=H_{\mathsf{S}}=\omega_{0}\sigma_{+}\sigma_{-}$. Hence from Eq. (20) of the main text, the short-time dynamics of the atom is obtained as

[TABLE]

from which it is seen that \varrho_{\mathrm{ee}}(\tau)=\varrho_{\mathrm{ee}}(0)+\tau^{2}\,\mathrm{Cov}_{\mathsf{B}_{0}}(\mathpzc{O}_{\mathsf{B}},\mathpzc{O}_{\mathsf{B}})\big{(}\varrho_{\mathrm{gg}}(0)-\varrho_{\mathrm{ee}}(0)\big{)}+O(\tau^{3}). Since the atom is assumed to be initially in the excited state $\varrho_{\mathrm{ee}}(0)=1$, thus

[TABLE]

which is in accordance with the short-time expansion of Eq. (21) of the main text. To compare Eq. (127) with the short-time expansion of the exact solution, we note that based on Eq. (9) of Ref. [22] of the main text and for the case of a single-mode bath $\langle\pm|\varrho_{\mathsf{S}}(\tau)|\mp\rangle=e^{-4f(\tau)}\langle\pm|\varrho_{\mathsf{S}}(0)|\mp\rangle$. Using the identity \varrho_{\mathrm{ee}}(\tau)=\frac{1}{2}\big{(}1-\langle+|\varrho_{\mathsf{S}}(\tau)|-\rangle-\langle-|\varrho_{\mathsf{S}}(\tau)|+\rangle\big{)}, which is obtained with a simple basis transformation using the relations $|+\rangle=(|\mathrm{g}\rangle+|\mathrm{e}\rangle)/\sqrt{2}$ and $|-\rangle=(|\mathrm{g}\rangle-|\mathrm{e}\rangle)/\sqrt{2}$, we conclude that $\varrho_{\mathrm{ee}}(\tau)=\frac{1}{2}(1+e^{-4f(\tau)})$, in which $f(\tau)=\tau^{2}\langle{\mathpzc{O}_{\mathsf{B}}^{2}}\rangle_{\mathsf{B}_{0}}/2$. Thus, the short-time behavior from the exact solution becomes $1-\tau^{2}\langle{\mathpzc{O}_{\mathsf{B}}^{2}}\rangle_{\mathsf{B}_{0}}+O(\tau^{3})$, which—considering $\langle{\mathpzc{O}_{\mathsf{B}}}\rangle_{\mathsf{B}_{0}}=0$—coincides with Eq. (127).

I.5 Comparison with the CG equation

In Fig. 7, we compare our data with the particular CG method delineated in Ref. Schaller-CG (Sec. III C 3 therein) with the MLL, ULL$2$, and other methods. We observe that here the MLL, ULL$2$, and TCL$2$ solutions outperform the CG and the Redfield solutions.

Appendix J Details of example III: Harmonic oscillator within a bath of oscillators

Here we provide the details related to example III introduced in Sec. VI B of the main text.

J.1 TCL$2$ equation

We recall the TCL$2$ equation (111). The interaction Hamiltonian for a harmonic oscillator interacting with a bath of similar oscillators is given by

[TABLE]

where $\hat{\bm{a}}(\tau)=\hat{a}e^{-i\omega_{0}\tau}$ and $\hat{\bm{b}}_{k}(\tau)=\hat{b}_{k}e^{-i\omega_{k}\tau}$ (see Ref. [39] of the main text for derivation). Now we can evaluate each term in the TCL$2$ equation. We have

[TABLE]

We assume that the bath is at zero temperature and simplify Eq. (129) further (using the relations $\mathrm{Tr}_{\mathsf{B}}[\hat{b}_{k}\hat{b}_{k^{\prime}}^{{\dagger}}\bm{\varrho}_{\mathsf{B}}(0)]=\delta_{kk^{\prime}}$, $\mathrm{Tr}_{\mathsf{B}}[\hat{b}_{k}^{{\dagger}}\hat{b}_{k^{\prime}}\bm{\varrho}_{\mathsf{B}}(0)]=0$, $\mathrm{Tr}_{\mathsf{B}}[\hat{b}_{k}\hat{b}_{k^{\prime}}\bm{\varrho}_{\mathsf{B}}(0)]=0$, and $\mathrm{Tr}_{\mathsf{B}}[\hat{b}_{k}^{{\dagger}}\hat{b}_{k^{\prime}}^{{\dagger}}\bm{\varrho}_{\mathsf{B}}(0)]=0$),

[TABLE]

Thus, the TCL$2$ equation reads as

[TABLE]

If we set $x=\tau-s$ we obtain

[TABLE]

where $\Delta_{k}=\omega_{0}-\omega_{k}$. Using the integral relation

[TABLE]

yields

[TABLE]

where $\gamma(\tau)=\sum_{k}g_{k}^{2}\sin(\Delta_{k}\tau)/\Delta_{k}$ and $\kappa(\tau)=\sum_{k}2g_{k}^{2}\sin^{2}(\Delta_{k}\tau/2)/\Delta_{k}$. After transforming back to the Schrödinger picture, we obtain the following TCL$2$ equation:

[TABLE]

As a remark, note that $\gamma(\tau)\approx G\tau$ for small $\tau$’s and $\lim_{\tau\rightarrow 0}\kappa(\tau)=0$, hence in this limit this TCL$2$ equation reduces to the MLL equation (23) of the main text.

J.2 CG equation

Here we use the CG method as per Ref. Majenz and obtain the following CG master equation for a harmonic oscillator within a bath of oscillators:

[TABLE]

where

[TABLE]

We have used the Ohmic spectral density $J(\omega)=(\omega/\pi)e^{-\omega/\omega_{c}}$, where $\omega_{c}$ is the cutoff frequency. The optimal value of the free parameter $\Delta t$ can be fixed by comparing with the exact solution.

J.3 ULL$2$ equation

We use the coupled ULL$2$ dynamical equations (106) and (108). The results of the numerical simulations of the the system dynamical equation have been shown in Fig. 4 of the main text. In Fig. 8 of this note we have depicted the population of the first excited state for example III for both ULL$2$ and exact solutions, when $\omega_{c}=10$ and the rest of the parameters have values equal to those used in Fig. 4 of the main text. Again, the ULL$2$ equation outperforms the other methods (see also the inset of Fig. 4 of the main text).