Succinct Description and Efficient Simulation of Non-Markovian Open Quantum Systems

Xiantao Li Department of Mathematics, Pennsylvania State University Chunhao Wang

Abstract

Non-Markovian open quantum systems represent the most general dynamics when the quantum system is coupled with a bath environment. The quantum dynamics arising from many important applications are non-Markovian. Although for special cases, such as Hamiltonian evolution and Lindblad evolution, quantum simulation algorithms have been extensively studied, efficient quantum simulations for the dynamics of non-Markovian open quantum systems remain underexplored. The most immediate obstacle for studying such systems is the lack of a universal succinct description of their dynamics. In this work, we fulfill the gap of studying such dynamics by 1) providing a succinct representation of the dynamics of non-Markovian open quantum systems with quantifiable error, and 2) developing an efficient quantum algorithm for simulating such dynamics with cost $\mathcal{O}(t\,\mathrm{polylog}(t/\epsilon))$ for evolution time $t$ and precision $\epsilon$ . Our derivation of the succinct representation is based on stochastic Schrödinger equations, which could lead to new alternatives to deal with open quantum systems as well.

1 Introduction

As the size of many modern-day electronic devices is continuously being reduced, quantum mechanical properties start to become dominant. Many novel designs have emerged, e.g., quantum wires, quantum dots, and molecular transistors, to take advantage of these properties. What is common in these applications is that the quantum properties are vital. However, these quantum dynamics do not evolve in isolation. Known as open quantum systems [BP02], they are always interacting with their environment, which due to the large dimension, can not be included explicitly in the computation. Another challenge comes from the fact that the continuous interactions can give rise to non-Markovian behavior, for which standard descriptions break down. The implication of non-Markovian dynamics to the quantum properties has been analyzed through many model examples [BLPV16, PZA⁺19, dVA17, SP16, SH04, AM21, MEL20]. More importantly, there are important experimental observations of non-Markovian dynamics [GTP⁺15, MALH⁺11], and such behavior has a strong impact on the electronic properties of molecular devices. Furthermore, it has been discovered that non-Markovianity can enhance quantum entanglement [TER⁺09].

Generally speaking, the dynamics of a quantum system interacting with a bath environment can be described by the von Neumann equation¹¹1The von Neumann equation is a generalization of the Schrödinger equation to the context of density matrices.

\displaystyle i\partial_{t}\rho=[H_{\text{tot}},\rho],\quad\rho(0)=\rho_{S}(0)\otimes\rho_{B},

(1)

where the Hamiltonian $H_{\text{tot}}$ that couples the system to a bath is expressed as,

H_{\text{tot}}=H_{S}\otimes I_{B}+I_{S}\otimes H_{B}+\lambda\sum_{\alpha=1}^{M}S_{\alpha}\otimes B_{\alpha}.

(2)

Here, $\lambda$ is typically referred to as the coupling parameter, and the integer $M$ indicates the number of interaction terms. To consider the problem in the context of quantum simulations, we let $2^{n}$ be the dimension of the system (S), i.e., $H_{S}\in\mathbb{C}^{2^{n}\times 2^{n}}$ is acting on $n$ qubits. A widely considered example is the spin Boson model, where $H_{S}$ and $S_{\alpha}$ are Pauli matrices. One notable example in this category, which is also relevant to quantum computers, is the dynamics of qubits coupled to a boson bath [WC13]. Other applications include the Anderson-Holstein model [CL17a] where $H_{S}$ consists of the kinetic and potential energy terms, and those from electronic transport problems, where the Hamiltonian is expressed in terms of molecular or atomic orbitals [GN05].

When $\lambda=0$ , the system is completely decoupled from the environment, and the dynamics of the system reduces to Hamiltonian evolution, which can also occur for finite $\lambda$ , but under more stringent assumptions on the operator $S_{\alpha}$ and the initial state [LCW98]. In recent decades, the problem of simulating Hamiltonian evolution has been extensively studied. A subset of notable works can be found in [Llo96, BCG14, BCC⁺15, BCK15, BCC⁺17, LC17, LC19, CGJ19]. This line of research has led to optimal Hamiltonian simulation algorithms [LC17, GSLW19].

In general, the continuous interactions with the bath conspire to a host of interesting quantum dynamics. Known as open quantum systems [BP02], such problems have led to a long-standing challenge in modeling the system dynamics without an explicit representation of the bath. Due to such modeling difficulty, the progress of developing fast simulation algorithms is much slower for open quantum systems. One exception is Markovian open quantum systems, which arises when $0<\lambda\ll 1$ , and the dynamics of the bath occurs at a much faster rate than the system. Intuitively, if an open quantum system is Markovian, the bath Hamiltonian is fast enough to “forget” the disturbance caused by the system-bath interaction, and therefore information only flows from system to bath with no transfer of information back to the system. By computing the infinitesimal generator, a complete description of the Markovian dynamics has been obtained [Lin76], which later is referred to as the Lindblad equation. In 2011, Kliesch, Barthel, Gogolin, Kastoryano, and Eisert [KBG⁺11] gave the first quantum algorithm for simulating Markovian open quantum systems with cost $\mathcal{O}(t^{2}/\epsilon)$ for evolution time $t$ and precision $\epsilon$ . In 2017, Childs and Li [CL17b] presented an algorithm that improves the cost to $\mathcal{O}(t^{1.5}/\sqrt{\epsilon})$ , and Cleve and Wang [CW17] further improved the algorithm by reducing the complexity to $\mathcal{O}(t\,\mathrm{polylog}(t/\epsilon))$ , which is nearly optimal. All these algorithms are designed for models with sparse Hamiltonian and jump operators. Another notable approach is by Schlimgen et al. [SHMS⁺21, SHMSS⁺22], where the coefficient matrices in the Kraus form associated with the Lindblad evolution are assumed to be given inputs. Overall, the difficulty with simulating Markovian open quantum systems lies in the observation that the advanced algorithmic techniques for Hamiltonian simulation cannot be directly applied to open quantum systems because of the presence of decoherence. In fact, Cleve and Wang [CW17] have shown that it is impossible to achieve linear dependence on evolution time by a direct reductionist approach.

In sharp contrast to quantum Markovian processes, there is no universal form for the quantum master equation for non-Markovian dynamics. In the regime when the rate of the bath dynamics is comparable to the rate of the system dynamics, the system becomes non-Markovian. Loosely speaking, the non-Markovian dynamics can be interpreted as the backflow of information for the environment to the open quantum system, and it can be more precisely characterized using various metrics [BLP09, VMP⁺11, LFS12]. For such systems, many approaches have been developed to describe such dynamics in a way that goes beyond the Lindblad framework [CK10, dVA17, DS97, GN99, IT05, KMCWM16, MT99, MCR16, MCR17, PMCKM19, PRRF⁺18, Sha04, SG03, SKK⁺18, Str96, SDG99, SES14, Tan06, Li21, Tan06]. The earliest approach dates back to the projection formalism of Nakajima and Zwanzig [Nak58, Zwa60]. At the level of density matrices, the non-Markovian property is reflected in a memory integral that involves the bath correlation function, which can be approximated by a linear combination of Lorentzian terms [RE14]. This representation enables an embedding procedure, where the memory integral can be replaced by the dynamics of additional density-matrices [MT99, Li21]. Although a unified framework is not yet present, such a procedure embeds the density matrix of the system in an extended system of equations involving additional density matrices, for which the dynamics is Markovian. The coefficients in the extended dynamics are connected to the spectral properties of the bath. Such a quantum master equation is often referred to as the generalized quantum master equation (GQME).

1.1 Main results

In this work, we 1) provide a succinct GQME representation for the dynamics of non-Markovian quantum systems, and 2) develop an efficient quantum algorithm for simulating the dynamics of open quantum systems based on this new representation. The derivation of the GQME is an important mathematical contribution that provides the groundwork for algorithmic development. It is worth noting that we only consider open quantum systems with coupling parameter $\lambda\ll 1$ . This is the scenario when we have a priori bound on the model error, whereas in the strong coupling regime, it is difficult to determine such an error in advance, and results can only be trusted with a leap of faith.

Mathematical contribution

In modeling open quantum systems, an important property to retain is the positivity. Toward this end, we choose to work with one that can be derived from an unraveling approach [BP02]. Namely, there exists an underlying stochastic Schrödinger equation (SSE) and the state-matrix is automatically positive semidefinite. In order to accurately incorporate the effect of the bath, let $K$ be the number of Lorentzian terms in representing the bath correlation function [RE14]. We consider an embedding of the system state $\rho_{S}$ into an $\mathcal{O}(n\log K)$ -qubit system (with dimension $K2^{n}$ ). Let $\Gamma$ denote the (unnormalized) state density matrix of this larger space. We derive from the SSE the following quantum master equation to describe the dynamics of $\Gamma$ :

\displaystyle\partial_{t}\Gamma=-i(H\Gamma-\Gamma H^{{\dagger}})+\sum_{k=1}^{K}V_{k}\Gamma V_{k}^{{\dagger}},

(3)

where $H$ (in general non-Hermitian) is defined in Eq. 41 which involves $H_{S}$ , $S_{\alpha}$ , and the bath correlation function (assumed to be part of the input). $V_{k}$ will be defined later in Eq. 39 which involves the bath correlation function. The initial state $\Gamma(0)$ , as deduced from the SSE, is given by,

\displaystyle\Gamma(0):=\outerproduct{0}{0}\otimes\rho_{S}(0)+\sum_{k=1}^{K}(\outerproduct{4k-2}{4k-2}+3\outerproduct{4k}{4k})\otimes\rho_{S}(0).

(4)

The system state is embedded in $\Gamma$ in the sense that the upper-left block of $\Gamma(t)$ — the unnormalized state of the larger system at time $t$ — approximates the exact system state $\rho_{S}(t)$ at time $t$ with error up to ${\mathcal{O}\left(\lambda^{3}\right)}$ (see Theorem 10). This error is referred to as the model error as it characterizes the accuracy of modelling the actual dynamics determined by Eq. 1 as a succinct GQME. Our model error of ${\mathcal{O}\left(\lambda^{3}\right)}$ improves the model error ${\mathcal{O}\left(\lambda^{3}/\alpha^{3}\right)}$ of the Lindblad equation with $\alpha\ll 1$ representing the time scale separation between the system and the bath, which (surprisingly) had not been precisely characterized until recently [CL17a].

One quick illustrative example of the GQME is a two-qubit model coupled to a common bosonic bath [WC13], where HEOM type of equations were derived. Specifically, the system Hamiltonian is written as $H_{s}=\frac{\omega_{0}}{2}(\sigma^{\text{I}}_{z}+\sigma^{\text{I\!I}}_{z})$ , where I and II label the two qubits and $\omega_{0}$ is the Zeeman energy. In addition, in the coupling term, $S=\sigma^{\text{I}}_{x}+\sigma^{\text{I\!I}}_{x}$ . Thus $M=1.$ Furthermore, the study in [WC13] considered the bath correlation function $C(t)=\frac{\lambda\gamma}{2}e^{-(\gamma+i\omega_{0})t}.$ In light of Eq. 5, we have that $K=1$ , $\theta_{1}=\sqrt{\frac{\lambda\gamma}{2}},$ $\ket{Q_{k}}=1,$ and $d_{1}=\omega_{0}+i\gamma.$ From the derivation in Eq. 36, we also have $V_{1}=\theta_{1}S.$ In this case, the matrix $\Gamma$ is a $5\times 5$ block matrix with total dimension being $20.$

Input model and simulation problem

The computational problem we consider is to simulate the dynamics generated by Eq. 3, which consists of an initial state preparation problem and a target state approximation problem. More formally, we formulate the problem as follows.

Problem 1 (Simulating non-Markovian open quantum systems).

Consider the dynamics defined by Eq. 3. Suppose we are given access to some efficient descriptions of the operators in Eq. 3. For any initial state $\rho_{S}(0)$ of the system, evolution time $t$ , and precision parameter $\epsilon$ , we need to

1.

prepare the initial state $\Gamma(0)$ as in Eq. 4, and
2.

produce a quantum state $\widetilde{\rho}_{S}(t)$ for the system so that the trace-distance between this state and the upper-left block of $\Gamma(t)$ , which is $\rho_{S}(t)$ , is at most $\epsilon$ , where $\Gamma(t)$ is the (unnormalized) state of evolving Eq. 3 for time $t$ with initial state $\Gamma(0)$ .

To solve this simulation problem, we need efficient descriptions of the operators in Eq. 3. The most straightforward input model is to assume that we are given these efficient descriptions directly. Such assumptions have been made in the literature of simulating Markovian open quantum systems [KBG⁺11, CL17b, CW17]. However, this straightforward input model is often not physically feasible: In many cases, we only have low-level information about the system, such as the system Hamiltonian and the system part of interaction Hamiltonians in Eq. 2, while high-level information such as descriptions of the operators in Eq. 3 is not readily available. With this practical consideration in mind, we work with a low-level input model, i.e., we only assume information that arises in Eq. 2, which is more convenient in real-world applications.

In particular, for the operators $H_{S}$ and $S_{\alpha}$ , we use a widely-used input model that has been recently introduced by Low and Chuang [LC19], and Chakraborty, Gilyén, and Jeffery [CGJ19] — the block-encodings of Hamiltonians. Roughly speaking, a block-encoding with normalizing factor $\alpha$ of a matrix $A$ is a unitary $U$ whose upper-left block is $A/\alpha$ . This input model is general enough to include almost all efficient representations that arise in physics applications, including linear combinations of tensor products of Paulis, sparse-access oracles, and local Hamiltonians. More specifically, if a matrix is $k$ -local, then it is $2^{k}$ -sparse. If a matrix $H$ is $d$ -sparse, then it can be approximated as a linear combination of unitaries with sum-of-coefficients $\mathcal{O}(d^{2}\norm{H}_{\max})$ , where $\norm{H}_{\max}$ is the largest entry of $H$ in absolute value. Moreover, if a matrix can be written as a linear combination of matrices with the sum-of-coefficients $s$ , then one can efficiently construct its block-encoding with normalizing constant $s$ . Note that the inverse directions of the above implications are in general not true. Therefore, our algorithm also works when a $H_{S}$ and $S_{\alpha}$ are provided in a less general input model, such as a sum of local operators, sparse-access oracles, or linear combinations of unitaries.

The BCF provides valuable information on how the bath influences the dynamics of the quantum system. In this work, we consider a typical representation of the bath correlation function (BCF), expressed as,

\displaystyle{C}_{\alpha,\beta}(t):=\sum_{k=1}^{K}\theta_{k}^{2}\innerproduct{\alpha}{Q_{k}}\innerproduct{Q_{k}}{\beta}\exp(-id_{k}^{*}t),

(5)

where $\ket{Q_{k}}\in\mathbb{C}^{M}$ , $\theta_{k}\in\mathbb{R}$ , and $d_{k}\in\mathbb{C}$ . Note that for all $t$ , ${C}_{\alpha,\beta}$ is an $M\times M$ matrix — a size that is tractable for classical computers. The treatment of the BCF is at the heart of modeling open quantum systems, and in practice, the specification of the BCF starts with the spectral density (SD) $J_{\alpha,\beta}(\omega)$ that depends on the bath spectrum and interaction strength. For bosonic environment, they are related as follows,

C_{\alpha,\beta}(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\big{(}\coth(\frac{\omega}{2k_{B}T})\cos(\omega t)-i\sin(\omega t)\big{)}J_{\alpha,\beta}(\omega)d\omega,

(6)

e.g., see [LRA⁺20], and [JZY08] for more general cases. Depending on the application, e.g., solvent, biomolecules, and nano materials, the spectral density can be adjusted accordingly. The pole expansion approach using Cauchy’s Residue Theorem has been applied to various types of SD functions [RE14]. In particular, both the functions $\coth(\frac{\omega}{2k_{B}T})$ and $J(\omega)$ can be treated this way. Such an expansion, after a truncation beyond a cut-off frequency $d_{\max}$ for the poles, gives rise to a finite sum of complex exponential terms [RE14]. The form of the coefficients in Eq. 5 is to ensure the positive definite property of the BCF. Namely, after the Fourier transform, the BCF has to be a semi positive definite function.

Our algorithm will take such approximation results as an input. We consider the poles $d_{k}$ within a cut-off frequency $d_{\max}$ , i.e., $\absolutevalue{d_{k}}\leq d_{\max}$ .

With the weak coupling condition, and the explicit representation Eq. 5 of the BCF, we will derive a generalized quantum master equation as in Eq. 3, with the extended quantum state $\Gamma(t)$ encoding the system density matrix $\rho_{S}(t).$ In addition, we demonstrate how the new Hamiltonian operator $H$ in Eq. 3, as well as the jump operators $V_{k}$ , can be obtained from the coefficients in the BCF Eq. 5.

In summary, in our quantum algorithm, we assume we are given the following quantities, as an efficient description of Eq. 3:

1.

a block-encoding $U_{H_{S}}$ of $H_{S}$ ;
2.

a block-encoding $U_{S_{\alpha}}$ of $S_{\alpha}$ for each $\alpha\in[M]$ ;
3.

the real numbers $\theta_{k}$ , vectors $\ket{Q_{k}}$ (all entries), and complex numbers $d_{k}$ for $k\in[K]$ as in Eq. 5, together with an upper bound $d_{\max}$ on $\absolutevalue{d_{k}}$ .

Algorithmic contribution

Our main algorithmic contribution is a quantum algorithm that solves 1. We informally state this result as follows.

Theorem 2 (Informal version of Theorem 18).

Suppose that we are given a block-encoding $U_{H_{S}}$ of $H_{S}$ , block-encodings $U_{S_{\alpha}}$ of $S_{\alpha}$ (for $\alpha\in[M]$ ), $\theta_{k}$ , $\ket{Q_{k}}$ (all entries), $\lambda$ , and $d_{k}$ for $k\in[K]$ as in Eq. 5. Then there exists a quantum algorithm that solves 1 using

\displaystyle{\mathcal{O}\left(t\,\mathrm{polylog}(t/\epsilon)\mathrm{poly}({\mu},M,K,\lambda,d_{\max})\right)},

(7)

queries to $U_{H_{S}}$ , and $U_{S_{a}}$ and additional 1- and 2-qubit gates, where $\mu$ is the normalizing factor of the block-encodings, $M$ is the number of interaction terms in Eq. 1, $K$ is the number of Lorentzian terms in Eq. 5, and $d_{\max}$ is an upper bound of $\absolutevalue{d_{k}}$ .

Our algorithm follows the high-level idea of [CW17], but we have generalized their techniques to work with block-encoded inputs. The building block of our algorithm is an implementation of completely positive maps given block-encodings of the Kraus operators (see Lemma 6 for more details). Suppose the normalizing factors of the block-encodings of the Kraus operators are $\alpha_{1},\ldots,\alpha_{m}$ , then our construction gives the success probability parameter $1/\sum_{j=1}^{m}\alpha_{j}^{2}$ , while a straightforward construction using Stinespring dilation yields a worse success probability parameter $1/(\sum_{j=1}^{m}\alpha_{j})^{2}$ , which does not permit the desired dependence on $t$ and $\epsilon$ .

Since the dynamics that Eq. 3 generates is completely positive, we consider an infinitesimal approximation map that approximates the first-order Taylor approximation of the dynamics for a small enough evolution time. From the low-level input model, we can construct the block-encodings of the Kraus operators of this infinitesimal approximation map, and it can be implemented using our building block Lemma 6 with high success probability parameter. We repeat this construction until the success probability parameter becomes a constant, and at that point, we have obtained a normalized version of $\Gamma(t)$ for $t$ proportional to a constant. Now, to extract the upper-left block of the resulting (normalized) density matrix, we use oblivious amplitude amplification for isometries ([CW17] and Lemma 17) to achieve this with an extra factor $\sqrt{K}$ , as the trace of $\Gamma(t)$ is upper bounded by $\mathcal{O}(K)$ for all $t$ (see Corollary 12).

For the problem of simulating Hamiltonian dynamics, high-order Taylor approximation yields simpler and faster quantum algorithms, e.g. [BCC⁺15]. However, for simulating open quantum systems, it is not known how to take advantage of higher-order approximations. This is because high powers of the superoperator defined in Eq. 3 are too complicated to keep track of its completely positive structure, which is the key to implementing such maps.

1.2 Summary of contributions

We highlight our contributions as follows:

1.

We provide a succinct representation of non-Markovian dynamics, where the density matrix from the GQME is consistent with that from the full quantum dynamics with provable ${\mathcal{O}\left(\lambda^{3}\right)}$ accuracy. A notable feature of our new succinct representation is that the positivity is guaranteed, which is the key requirement for designing quantum simulation algorithms.
2.

We develop a quantum algorithm based on this representation. The cost of our algorithm scales linearly in $t$ , poly-logarithmically in $\epsilon$ , and polynomially in $M$ and $K$ . To the best of our knowledge, this algorithm is the first to achieve linear dependence on $t$ and poly-logarithmic dependence on $\epsilon$ for simulating non-Markovian open quantum systems. In addition, our algorithm works with low-level input models, which are readily available in many real-world applications.
3.

Other technical contributions: We have shown that the GQMEs can be unraveled (see [BP02]) into stochastic Schrödinger equations, which provides another potential alternative to obtain the density matrix as a statistical quantity. In addition, we prove that the extended density matrix from the GQME is bounded over the time scale ${\mathcal{O}\left(\lambda^{-1}\right)}.$

1.3 Related work

In the context of modeling open quantum systems [BP02], the hierarchical equations of motion (HEOM) approach [Tan06, Tan20] also yields an extended dynamics for the density matrix, but without using the weak coupling assumption. Rather, the equations are truncated based on numerical observations. In principle, the quantum algorithms in this work can also be applied to those GQMEs from the HEOM approach, provided that the approach can be proved to produce completely positive maps. Meanwhile, since no error bound is available, it is difficult to prescribe the level of truncation in advance, which also makes it difficult to estimate the computational complexity.

Another interesting development is to embody the memory effect using a local form of the GQME, where the generator consists of multiple Lindblad operators with time-dependent coefficients. In fact, it has been proved [CK10] that any nonlocal form of the non-Markovian dynamics can be rewritten in a local form with time-dependent generators. Due to the fact that the proof is non-constructive, the implementations of this type of Lindblad operators are empirical. Sweke, Sanz, Sinayskiy, Petruccione, and Solano [SSS⁺16] considered such time-local quantum master equations, and constructed quantum algorithms. Their algorithms rely on Trotter splittings, by decomposing the entire generator into local operators. Since it is not yet clear how the GQMEs from the current approach, or those from the HEOM approach, can be expressed in time-local forms, a direct comparison of the computational complexity is not yet available.

1.4 Open questions

Modeling open quantum systems outside the weak coupling regime is still an outstanding challenge. The HEOM approach [Tan20] relies on a frequency cut-off to achieve a Markovian embedding. Quantifying the error associated with such an approximation, and ensuring the positivity are two of the remaining issues.

Another interesting scenario is when the open quantum system is subject to an external potential. Quantum optimal control is one important example. Deriving a GQME in the presence of a time-dependent external field while still maintaining the control properties is still an open problem to the best of our knowledge.

The cost of our quantum algorithm is $\mathcal{O}(t\,\mathrm{polylog}(t/\epsilon))$ . Is there a faster quantum algorithm that achieves an additive cost, i.e., $\mathcal{O}(t+\mathrm{polylog}(1/\epsilon))$ ? This additive cost is the lower bound for Hamiltonian simulation [BACS07, BCK15, BCC⁺17], and it is hence a lower bound for the non-Markovian simulation problem. The optimal Hamiltonian simulation was achieved by quantum signal processing due to Low and Chuang [LC17], which has been generalized to quantum singular value transformation by Gilyén, Su, Low, and Wiebe [GSLW19]. Unfortunately, these techniques do not immediately generalize to open quantum systems, as it is not clear what the correspondence of singular values and eigenvalues should be for superoperators.

2 Preliminaries

2.1 Notation

In this paper, we use the ket-notation to denote a vector only when it is normalized. For a vector $v$ , we use $\norm{v}$ to denote its Euclidean norm. For a square matrix $M$ , we use $\norm{M}$ to denote its spectral norm and use $\norm{M}_{1}$ to denote its trace norm, i.e., $\norm{M}_{1}=\tr(\sqrt{M^{{\dagger}}M})$ . The identity operator acting on a Hilbert space of dimension $N$ is denoted by $I_{N}$ , e.g., $I_{2^{n}}$ is the identity operator acting on $n$ qubits. When the context is clear, we drop the subscript and simply use $I$ . We use calligraphic fonts, such as $\mathcal{K}$ , $\mathcal{L}$ , and $\mathcal{M}$ to denote superoperators, which maps matrices to matrices. We consider superoperators of the form

\displaystyle\mathcal{M}:\mathbb{C}^{N\times N}\rightarrow\mathbb{C}^{M\times M},

(8)

and use $\mathrm{T}(\mathbb{C}^{N},\mathbb{C}^{M})$ to denote the set of all such superoperators. In particular, we use $\mathcal{I}$ to denote the identity map, which maps every matrix to itself. Whenever necessary, we use the subscript $N$ of $\mathcal{I}_{N}\in\mathrm{T}(\mathbb{C}^{N},\mathbb{C}^{N})$ to highlight the dimension of matrices it acts on. For example, $\mathcal{I}_{2^{n}}$ is acting on $n$ -qubit operators. The induced trace norm of a superoperator $\mathcal{M}\in\mathrm{T}(\mathbb{C}^{N},\mathbb{C}^{M})$ , denoted by $\norm{\mathcal{M}}_{1}$ , is defined as

\displaystyle\norm{\mathcal{M}}_{1}=\max\{\norm{\mathcal{M}(A)}_{1}:A\in\mathbb{C}^{N\times N},\norm{A}_{1}\leq 1\}.

(9)

The diamond norm of a superoperator $\mathcal{M}\in\mathrm{T}(\mathbb{C}^{N},\mathbb{C}^{M})$ , denoted by $\norm{\mathcal{M}}_{\diamond}$ is defined as

\displaystyle\norm{\mathcal{M}}_{\diamond}:=\norm{\mathcal{M}\otimes\mathcal{I}_{N}}_{1}.

(10)

For a positive integer $m$ , we use $[m]$ to denote the set $\{1,2,\ldots,m\}$ .

2.2 The block-encoding method

We use the notion of block-encoding as an efficient description of input operators. To have access to an operator $A$ , we assume we have access to a unitary $U$ whose upper-left block encodes $A$ in the sense that

\displaystyle U=\begin{pmatrix}A/\alpha&\cdot\\ \cdot&\cdot\\ \end{pmatrix},

(11)

which implies that $A=\alpha(\bra{0}\otimes I)U\ket{0}\otimes I)$ . More precisely, we have the following definition.

Definition 3.

Let $A$ be an $n$ -qubit operator. For a positive real number $\alpha>0$ and natural number $m$ , we say that an $(n+m)$ -qubit unitary $U$ is an $(\alpha,m,\epsilon)$ -block-encoding of $A$ if

\displaystyle\norm{A-\alpha(\bra{0}\otimes I_{2^{n}})U\ket{0}\otimes I_{2^{n}})}\leq\epsilon.

(12)

The following lemma shows how to construct a block-encoding for sparse matrices.

Lemma 4 ([GSLW19, Lemma 48]).

Let $A\in\mathbb{C}^{2^{n}\times 2^{n}}$ be an $n$ -qubit operator with at most $s$ nonzero entries in each row and column. Suppose $A$ is specified by the following sparse-access oracles:

	$\displaystyle O_{A}:$	$\displaystyle\ket{i}\ket{j}\ket{0}\mapsto\ket{i}\ket{j}\ket{A(i,j)},\text{ and}$		(13)
	$\displaystyle O_{S}:$	$\displaystyle\ket{i}\ket{k}\mapsto\ket{i}\ket{r_{i,k}},$		(14)

where $r_{i,k}$ is the $k$ -th nonzero entry of the $i$ -th row of $A$ . Suppose $\absolutevalue{A_{i,j}}\leq 1$ for $i\in[m]$ and $j\in[n]$ . Then for all $\epsilon\in(0,1)$ , an $(s,n+3,\epsilon)$ -block-encoding of $A$ can be implemented using $\mathcal{O}(1)$ queries to $O_{A}$ and $O_{S}$ , along with $\mathcal{O}(n+\mathrm{polylog}(1/\epsilon))$ 1- and 2-qubit gates. Moreover, if $A_{i,j}\in\{0,1\}$ for all $i\in[m]$ and $j\in[n]$ , the block-encoding of $A$ can be implemented precisely, i.e., $\epsilon=0$ .²²2The case when $A_{i,j}\in\{0,1\}$ was not explicitly stated in [GSLW19, Lemma 48]; however, the conclusion is not hard to obtain as a special case of their proof.

A linear combination of block-encodings can be constructed using [GSLW19, Lemma 52]. Here, we slightly generalize their construction to achieve better performance when the normalizing factors of each block-encoding are different.

Lemma 5.

Suppose $A:=\sum_{j=1}^{m}y_{j}A_{j}\in\mathbb{C}^{2^{n}\times 2^{n}}$ , where $A_{j}\in\mathbb{C}^{2^{n}\times 2^{n}}$ and $y_{j}>0$ for all $j\in\{1,\ldots m\}$ . Let $U_{j}$ be an $(\alpha_{j},a,\epsilon)$ -block-encoding of $A_{j}$ , and $B$ be a unitary acting on $b$ qubits (with $m\leq 2^{b}-1$ ) such that $B\ket{0}=\sum_{j=0}^{2^{b}-1}\sqrt{\alpha_{j}y_{j}/s}\ket{j}$ , where $s=\sum_{j=1}^{m}y_{j}\alpha_{j}$ . Then a $(\sum_{j}y_{j}\alpha_{j},a+b,\sum_{j}y_{j}\alpha_{j}\epsilon)$ -block-encoding of $\sum_{j=1}^{m}y_{j}A_{j}$ can be implemented with a single use of $\sum_{j=0}^{m-1}\outerproduct{j}{j}\otimes U_{j}+((I-\sum_{j=0}^{m-1}\outerproduct{j}{j})\otimes I_{\mathbb{C}^{2^{a}}}\otimes I_{\mathbb{C}^{2^{n}}})$ plus twice the cost for implementing $B$ .

Proof.

The proof is similar to that of [GSLW19, Lemma 52]. The difference is that, instead of preparing the state $\sum_{j=1}^{m}\frac{\sqrt{y_{j}}}{\sqrt{\sum_{j=1}^{m}y_{j}}}\ket{j}$ , we use the state preparation gate $B$ here. First note that

\displaystyle\norm{(\bra{0}^{\otimes a}\otimes I_{2^{n}})U_{j}(\ket{0}^{\otimes a}\otimes I_{2^{n}})-A_{j}/\alpha_{j}}\leq\epsilon.

(15)

Let $W=\sum_{j=0}^{m-1}\outerproduct{j}{j}\otimes U_{j}+((I-\sum_{j=0}^{m-1}\outerproduct{j}{j})\otimes I_{2^{a}}\otimes I_{2^{n}})$ and define $\widetilde{W}=(B^{{\dagger}}\otimes I_{2^{a}}\otimes I_{2^{n}})W(B\otimes I_{2^{a}}\otimes I_{2^{n}})$ . We have

	$\displaystyle\norm{\sum_{j=1}^{m}y_{j}A_{j}-s\left(\bra{0}^{\otimes b}\otimes\bra{0}^{\otimes a}\otimes I_{2^{n}}\right)\widetilde{W}\left(\ket{0}^{\otimes b}\otimes\ket{0}^{\otimes a}\otimes I_{2^{n}}\right)}$		(16)
	$\displaystyle=\norm{\sum_{j=1}^{m}y_{j}A_{j}-\sum_{j=1}^{m}\alpha_{j}y_{j}\left(\bra{0}^{\otimes a}\otimes I_{2^{n}})U_{j}(\ket{0}^{\otimes a}\otimes I_{2^{n}})\right)}$		(17)
	$\displaystyle\leq\sum_{j=1}^{m}y_{j}\alpha_{j}\norm{A_{j}/\alpha_{j}-\left(\bra{0}^{\otimes a}\otimes I_{2^{n}})U_{j}(\ket{0}^{\otimes a}\otimes I_{2^{n}})\right)}$		(18)
	$\displaystyle\leq\sum_{j=1}^{m}y_{j}\alpha_{j}\epsilon.$		(19)

∎

2.3 Technical tool for implementing completely positive maps

In this subsection, we provide the technical primitives for developing a simulation algorithm. The following lemma generalizes the technique of linear combination of unitaries for completely positive maps [CW17] to the context of block-encoding. This tool might be of independent interest as well.

Lemma 6.

Let $A_{1},\ldots,A_{m}\in\mathbb{C}^{2^{n}}$ be the Kraus operators of a completely positive map $\mathcal{M}$ [Kra71],

\mathcal{M}(\rho)=\sum_{j=1}^{m}A_{j}\rho A_{j}^{\dagger}.

Let $U_{1},\ldots,U_{m}\in\mathbb{C}^{2^{n+n^{\prime}}}$ be their corresponding $(s_{j},n^{\prime},\epsilon)$ -block-encodings, i.e.,

\displaystyle\norm{A_{j}-s_{j}(\bra{0}\otimes I)U_{j}\ket{0}\otimes I)}\leq\epsilon,\quad\text{ for all $1\leq j\leq m$}.

(20)

Let $\ket{\mu}:=\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j=1}^{m}s_{j}\ket{j}$ . Then $(\sum_{j=1}^{m}\outerproduct{j}{j}\otimes U_{j})\ket{\mu}\ket{0}\otimes I$ implements this completely positive map $\mathcal{M}$ in the sense that

\displaystyle\norm{I\otimes\bra{0}\otimes I\left(\sum_{j=1}^{m}\outerproduct{j}{j}\otimes U_{j}\right)\ket{\mu}\ket{0}\ket{\psi}-\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j=1}^{m}\ket{j}A_{j}\ket{\psi}}\leq\frac{m\epsilon}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}

(21)

for all $\ket{\psi}$ .

Proof.

It is easy to verify that

\displaystyle\left(\sum_{j=1}^{m}\outerproduct{j}{j}\otimes U_{j}\right)\ket{\mu}\ket{0}\ket{\psi}=\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j}s_{j}\ket{j}U_{j}(\ket{0}\ket{\psi}).

(22)

With a direct substitution into Eq. 21, one arrives at,

	$\displaystyle\norm{(I\otimes\bra{0}\otimes I)\left(\sum_{j=1}^{m}\outerproduct{j}{j}\otimes U_{j}\right)\ket{\mu}\ket{0}\ket{\psi}-\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j}\ket{j}A_{j}\ket{\psi}}$		(23)
	$\displaystyle=\norm{(I\otimes\bra{0}\otimes I)\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j}s_{j}\ket{j}U_{j}(\ket{0}\ket{\psi})-\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\sum_{j}\ket{j}A_{j}\ket{\psi}}$		(24)
	$\displaystyle=\frac{1}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}\norm{\sum_{j}s_{j}\ket{j}(\bra{0}\otimes I)U_{j}(\ket{0}\otimes I)\ket{\psi}-\sum_{j}\ket{j}A_{j}\ket{\psi}}$		(25)
	$\displaystyle\leq\frac{m\epsilon}{\sqrt{\sum_{j=1}^{m}s_{j}^{2}}}.$		(26)

∎

3 Non-Markovian Quantum Master Equation

In this section, we present the derivation of the GQME [Li21]. Non-Markovian dynamics has been extensively studied in the context of open quantum systems [BP02]. The starting point for considering an open quantum system is a quantum dynamics that couples a quantum system and a bath environment,

i\partial_{t}\rho=[H_{\text{tot}},\rho],\quad\rho(0)=\rho_{S}(0)\otimes\rho_{B},

(27)

where the coupled Hamiltonian is given by,

H_{\text{tot}}=H_{S}\otimes I_{B}+I_{S}\otimes H_{B}+\lambda\sum_{\alpha=1}^{M}S_{\alpha}\otimes B_{\alpha}.

(28)

Here $H_{S}\in\mathbb{C}^{2^{n}\times 2^{n}}$ is acting on $n$ qubits, and $M$ refers to the number of interaction terms.

We follow the standard setup by choosing $\rho_{B}$ according to a statistical ensemble, e.g.,

\rho_{B}=\frac{\exp(-\frac{H_{B}}{k_{B}T})}{Z},\quad Z:=\tr\left(\exp(-\frac{H_{B}}{k_{B}T})\right),

(29)

with $k_{B}$ and $T$ being respectively the Boltzmann constant and the temperature. As a result, $[H_{B},\rho_{B}]=0.$ Without loss of generality, we can assume that $\text{tr}(\rho_{B}B_{\alpha})=0,\;\alpha=1,2,\cdots,M,$ which can be ensured by properly shifting of the operators [GN99]. This helps to eliminate $\mathcal{O}(\lambda)$ terms in the asymptotic expansion [GN99].

Of primary interest in the theory of open quantum systems is the density matrix of the system, $\rho_{S}(t)$ , which in principle can be obtained with a partial trace: $\rho_{S}(t)=\text{tr}_{B}\left(\rho(t)\right).$ To arrive at a quantum master equation that embodies non-Markovian properties, we start with the non-Markovian stochastic Schrödinger equation (NMSSE), which was derived from the wave function representation of Eq. 27 in [GN99], and later revisited in [BD12]. In NMSSE, a stochastic realization of the quantum state follows a stochastic differential equation,

i\partial_{t}\psi=\hat{H}_{S}\psi-i\lambda^{2}\sum_{\alpha,\beta=1}^{M}\int_{0}^{t}C_{\alpha,\beta}(\tau)S^{\dagger}_{\alpha}e^{-i\hat{H}_{S}\tau}S_{\beta}\psi(t-\tau)d\tau+\lambda\sum_{\beta=1}^{M}\eta_{\beta}(t)S_{\beta}\psi(t).

(30)

Here $i=\sqrt{-1}$ . The matrix $C(t):\mathbb{R}\to\mathbb{C}^{M\times M}$ with elements $C_{\alpha,\beta}(t)$ corresponding to the correlation among the bath correlation $\{B_{\alpha}\}_{1\leq\alpha\leq M}$ .

Each noise term $\eta_{\alpha}(t):\mathbb{R}\to\mathbb{C}$ is Gaussian with mean zero and correlation given by [BD12],

\mathbb{E}[\eta_{\alpha}^{*}(t)\eta_{\beta}(t^{\prime})]=C_{\alpha,\beta}(t-t^{\prime}),\quad 1\leq\alpha,\beta\leq M.

(31)

The stationarity of the process also implied that $C(t)=C(-t)^{\dagger}$ . Thus, it suffices to consider the correlation function for $t\geq 0.$ This relation between a dissipation kernel and the time correlation of the noise is well-known in non-equilibrium statistical physics, and it is often labeled as the second fluctuation-dissipation theorem [Kub66].

3.1 The generalized quantum master equation (GQME)

With a scale separation assumption, the NMSSE can be reduced to the Lindblad equation [CL17a, LBW01]. In this regime, the bath correlation behaves like a delta function [GN99],

C(t-t^{\prime})\approx\sum_{j=1}^{M}\theta_{j}^{2}\outerproduct{R_{j}}{R_{j}}\delta(t-t^{\prime}),

(32)

which simplifies the memory integrals to local terms. By defining operators $L_{j}$ ,

{L_{j}}=\sum_{\beta=1}^{M}\theta_{j}\innerproduct{R_{j}}{\beta}S_{\beta},

(33)

the NMSSE then implies the following Lindblad equation [BD12],

i\partial_{t}\rho_{S}=[H_{S},\rho_{S}]-\lambda^{2}{\sum_{j=1}^{M}\Big{(}L^{\dagger}_{j}L_{j}\rho_{S}+\rho_{S}L^{\dagger}_{j}L_{j}-2L_{j}\rho_{S}L^{\dagger}_{j}\Big{)}.}

(34)

However, in the non-Markovian regime, the assumption Eq. 32 breaks down, i.e., the correlation length is finite, and a different approach is needed to derive a quantum master equation that governs the dynamics of $\rho_{S}$ . In this paper, we follow the general unraveling approach [BP02]. In Appendix A we show that the memory effect can be embedded in an extended stochastic system by introducing auxiliary wave functions, $\chi^{\text{I}}$ to $\chi^{\text{I\!V}}.$ The dynamics can be summarized in the following form,

\left\{\begin{array}[]{l}i\partial_{t}\psi=\displaystyle\hat{H}_{S}\psi-i\lambda\sum_{k=1}^{K}T_{k}^{\dagger}\chi^{\text{I}}_{k}-i\lambda\sum_{k}T_{k}\chi^{\text{I\!I}}_{k},\\ \begin{array}[]{ll}i\partial_{t}\chi^{\text{I}}_{k}=&(\hat{H}_{S}+d^{*}_{k})\chi^{\text{I}}_{k}+i\lambda T_{k}\psi(t),\\ i\partial_{t}\chi^{\text{I\!I}}_{k}=&(\hat{H}_{S}-d_{k})\chi^{\text{I\!I}}_{k}-i\lambda T_{k}^{\dagger}\chi^{\text{I\!I\!I}}_{k}-i\lambda T_{k}\chi^{\text{I\!V}}_{k}+i\gamma_{k}\psi(t)\dot{w}_{k}(t),\\ i\partial_{t}\chi^{\text{I\!I\!I}}_{k}=&i\lambda T_{k}\chi_{k}^{\text{I\!I}}+(\hat{H}_{S}-d_{k}-d_{k}^{*})\chi^{\text{I\!I\!I}}_{k}+i\gamma_{k}\chi_{k}^{\text{I}}\dot{w}_{k}(t),\\ i\partial_{t}\chi^{\text{I\!V}}_{k}=&i\lambda T_{k}^{\dagger}\chi_{k}^{\text{I\!I}}+(\hat{H}_{S}-2d_{k})\chi^{\text{I\!V}}_{k}+2i\gamma_{k}\chi_{k}^{\text{I\!I}}\dot{w}_{k}(t),\\ \end{array}\end{array}\right.

(35)

for $k=1,2,\ldots,K.$ Due to the presence of the auxiliary wave functions, the dynamics of $\psi$ is non-Markovian, and the additional equations induce a memory effect that imitates the non-Markovian behavior. We have dropped $\mathcal{O}(\lambda)$ terms in the last two equations, which can be justified as follows, by a substitution into the third equation, one can see that this truncation will contribute to an $\mathcal{O}(\lambda^{2})$ error to the dynamics of $\chi^{\text{I\!I}},$ which, after another substitution, leads to an $\mathcal{O}(\lambda^{3})$ perturbation in the first equation in Eq. 35. Namely, the accuracy is the same as the NMSSE Eq. 30. The operators $T_{k}$ are defined in terms of the operators $S_{\beta}$ in Eq. 2 and the coefficients in the BCF Eq. 5 as follows,

T_{k}=\sum_{\beta=1}^{M}\theta_{k}\innerproduct{Q_{k}}{\beta}S_{\beta},

(36)

From the definitions of the auxiliary functions, we can deduce their initial conditions,

\chi_{k}^{\text{I}}(0)=0,\quad\chi_{k}^{\text{I\!I}}(0)=i\zeta_{k}(0)\psi(0),\quad\chi_{k}^{\text{I\!I\!I}}(0)=0,\quad\chi_{k}^{\text{I\!V}}(0)=-\zeta_{k}(0)^{2}\psi(0),

with $\zeta_{k}$ ’s being independent Gaussian random variables of zero mean and unit variance.

To derive the corresponding GQME of the stochastic model in Eq. 116, we first write it as a system of SDEs,

i\partial_{t}\Psi=H\Psi+\sum_{k=1}^{K}V_{k}\Psi\dot{w}_{k}.

(37)

Here the function $\Psi$ includes the wave function $\psi$ and all the auxiliary wave functions $\{\chi_{k}^{\text{I}},\chi_{k}^{\text{I\!I}},\chi_{k}^{\text{I\!I\!I}},\chi_{k}^{\text{I\!V}}\}_{k=1}^{K}$ . One can arrange the wave functions and the operator $V_{1}$ as follows,

\Psi=\left(\begin{array}[]{c}\psi\\ \chi_{1}^{\text{I}}\\ \chi_{1}^{\text{I\!I}}\\ \chi_{1}^{\text{I\!I\!I}}\\ \chi_{1}^{\text{I\!V}}\\ \chi_{2}^{\text{I}}\\ \chi_{2}^{\text{I\!I}}\\ \chi_{2}^{\text{I\!I\!I}}\\ \chi_{2}^{\text{I\!V}}\\ \vdots\end{array}\right),\;V_{1}=\left(\begin{array}[]{ccccccc}0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&\dots\\ \sqrt{2\nu_{1}}&0&0&0&0&0&\dots\\ 0&\sqrt{2\nu_{1}}&0&0&0&0&\dots\\ 0&0&2\sqrt{2\nu_{1}}&0&0&0&\dots\\ 0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&\dots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)\otimes I_{S}.\;

(38)

Similarly,

V_{2}=\left(\begin{array}[]{ccccccccc}0&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&0&0&0&\dots\\ \sqrt{2\nu_{2}}&0&0&0&0&0&0&0&\dots\\ 0&0&0&0&0&\sqrt{2\nu_{2}}&0&0&\dots\\ 0&0&0&0&0&0&2\sqrt{2\nu_{2}}&0&\dots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)\otimes I_{S}.

Clearly, $V_{k}$ ’s are sparse and low rank. Consider the standard basis in $\mathbb{R}^{4K+1}$ , here abbreviated simply into $\ket{0},\ket{1},\ldots,\ket{4K}$ . A careful inspection reveals the following closed-form expressions,

\displaystyle V_{k}=\sqrt{2\nu_{k}}\Big{(}\outerproduct{{4k-2}}{0}+\outerproduct{{4k-1}}{{4k-3}}+2\outerproduct{{4k}}{{4k-2}}\Big{)}\otimes I_{S},

(39)

for $k=1,2,\ldots,K.$ Here $\nu_{k}=\text{Im}d_{k},$ and they are nonnegative.

The operator $H$ in Eq. 37 can also be written in a block form,

\left[\begin{array}[]{cccccccccc}H_{s}&-i\lambda T_{1}^{\dagger}&-i\lambda T_{1}&0&0&-i\lambda T_{2}^{\dagger}&-i\lambda T_{2}&0&0&\cdots\\ i\lambda T_{1}&H_{s}\!+\!d_{1}^{*}&0&0&0&0&0&0&0&\cdots\\ 0&0&H_{s}\!-\!d_{1}&-i\lambda T_{1}^{\dagger}&-i\lambda T_{1}&0&0&0&0&\cdots\\ 0&0&i\lambda T_{1}&H_{s}\!-\!2i\nu_{1}&0&0&0&0&0&\cdots\\ 0&0&i\lambda T_{1}^{\dagger}&0&H_{s}\!-\!2d_{1}&0&0&0&0&\cdots\\ i\lambda T_{2}&0&0&0&0&H_{s}\!+\!d_{2}^{*}&0&0&0&\cdots\\ 0&0&0&0&0&0&H_{s}\!-\!d_{2}&-i\lambda T_{2}^{\dagger}&-i\lambda T_{2}&\cdots\\ 0&0&0&0&0&0&i\lambda T_{2}&H_{s}\!-\!2i\nu_{2}&0&\cdots\\ 0&0&0&0&0&0&i\lambda T_{2}&0&H_{s}\!-\!2d_{2}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right].

(40)

The extended stochastic dynamics Eq. 35 introduced an auxiliary space that mimics the effect of the quantum bath. Recall $I_{S}\in\mathbb{C}^{2^{n}\times 2^{n}}$ is the identity operator. Similarly, we let $I_{A}\in\mathbb{C}^{(4K+1)\times(4K+1)}$ be the identify operator in an auxiliary space labelled by $A$ . Then the Hamiltonian in Eq. 40 can be expressed in terms of tensor products,

H=I_{A}\otimes H_{S}+H_{A}\otimes I_{S}+i\lambda\sum_{k=1}^{K}D_{k}\otimes T_{k}+i\lambda\sum_{k=1}^{K}E_{k}\otimes T_{k}^{\dagger}.

(41)

Here $H_{A}\in\mathbb{C}^{(4K+1)\times(4K+1)}$ is a diagonal matrix:

H_{A}=\text{diag}\left\{0,d_{1}^{*},-d_{1},d_{1}^{*}-d_{1},-2d_{1},d_{2}^{*},-d_{2},d_{2}^{*}-d_{2},-2d_{2},\ldots\right\}.

In addition, the matrices $D_{k},E_{k}\in\mathbb{C}^{(4K+1)\times(4K+1)}$ are given by

	$\displaystyle D_{k}=$	$\displaystyle\outerproduct{{4k-3}}{0}-\outerproduct{0}{{4k-2}}+\outerproduct{{4k-1}}{{4k-2}}-\outerproduct{{4k-2}}{{4k-1}},$		(42)
	$\displaystyle E_{k}=$	$\displaystyle-\outerproduct{0}{{4k-3}}+\outerproduct{{4k}}{{4k-2}}-\outerproduct{{4k-2}}{{4k}}.$		(42)

Assuming that the dimension of the original wave function is $n$ , i.e., $\ket{\psi}\in\mathbb{C}^{n}$ , the dimension of $\Psi$ is $(4K+1)n.$ Hence, the dimension of $\Gamma$ is $[(4K+1)n]\times[(4K+1)n]$ .

Within the framework of quantum unravelling [BP02], the density matrix associated with the combined wave functions $\Psi$ in Eq. 38 is defined as the entry-wise expectation,

\Gamma_{\alpha,\beta}(t)=\mathbb{E}[\Psi_{\alpha}(t)\Psi_{\beta}(t)^{*}].

(43)

An application of the Itô’s formula [Pav14] yields the following close-form quantum master equation,

\partial_{t}\Gamma=\mathcal{K}(\Gamma):=-i(H\Gamma-\Gamma H^{\dagger})+\sum_{k=1}^{K}V_{k}\Gamma V_{k}^{\dagger}.

(44)

The noise has been averaged out by the expectation. In fact, it has been shown in [KP13] that for any linear SDEs, the first and second moments satisfy close-form equations. Intuitively, the Lindblad description breaks down when the dynamics of $\rho_{S}(t)$ exhibits strong memory effect, i.e., it depends on the past history of $\rho_{S}(t)$ , which can not be captured by Eq. 34. On the other hand, the GQME Eq. 44 embodies the history dependence by embedding the dynamics of $\rho_{S}(t)$ in a larger system. Conceptually, if one solves other components of $\Gamma(t)$ and substitutes those solutions to the first block, one would get the memory dependence on $\rho_{S}(t)$ .

We can deduce the initial condition of $\Gamma$ from the definitions of the auxiliary wave functions. Since $\zeta_{k}$ is Gaussian with mean zero and variance 1, we have

\mathbb{E}[\innerproduct{\chi_{k}^{\text{I}}}{\chi_{k}^{\text{I}}}]=0,\;\mathbb{E}[\innerproduct{\chi_{k}^{\text{I\!I}}}{\chi_{k}^{\text{I\!I}}}]=\rho_{s}(0),\;\mathbb{E}[\innerproduct{\chi_{k}^{\text{I\!I\!I}}}{\chi_{k}^{\text{I\!I\!I}}}]=0,\;\mathbb{E}[\innerproduct{\chi_{k}^{\text{I\!V}}}{\chi_{k}^{\text{I\!V}}}]=3\rho_{s}(0).

All the cross-correlations are zero. Therefore, the initial density matrix $\Gamma$ is a block-diagonal matrix,

\Gamma(0)=\text{diag}\big{\{}\rho_{S}(0),0,\rho_{S}(0),0,3\rho_{S}(0),0,\rho_{S}(0),0,3\rho_{S}(0),\cdots\big{\}}.

(45)

We notice that the trace of $\Gamma(0)$ is $4K+1$ .

3.2 Properties of the GQME

We first provide some basic estimates on the extended density matrix $\Gamma$ . We begin by writing the Hamiltonian in Eq. 40 as,

H=H_{0}+\lambda H_{1}.

(46)

Here we made the observation that $H_{0}$ is block diagonal. On the other hand, $H_{1}$ only contains nonzero blocks in the first row and the first column. To conveniently refer to the block entries of the density matrix $\Gamma$ , we write it in the following block form,

\Gamma=\left[\begin{array}[]{cccc}\Gamma_{0,0}&\Gamma_{0,1}&\Gamma_{0,2}&\cdots\\ \Gamma_{1,0}&\Gamma_{1,1}&\Gamma_{1,2}&\cdots\\ \vdots&\vdots&\vdots&\ddots\\ \Gamma_{4K,0}&\Gamma_{4K,1}&\Gamma_{4K,2}&\cdots\end{array}\right].

(47)

In particular, $\rho_{S}$ is embedded into $\Gamma$ as the first block: $\rho_{S}=\Gamma_{0,0}.$

For such block matrices, we will use the following induced norm,

\norm{\Gamma}_{\infty}:=\max_{0\leq j\leq 4K}\sum_{0\leq i\leq 4K}\norm{\Gamma_{i,j}}.

(48)

Namely, for each entry, we use the spectral norm. But among the blocks, we use the $\infty$ -norm. One can verify that this norm still has the submultiplicative property. We choose this norm merely because we will estimate the bound of each block, and the formula in Eq. 48 can easily connect such estimates to the bound of the entire matrix. In principle, since matrix norms are continuous with respect to the entries, one can also use other norms among the blocks.

Lemma 7.

The exponential operator $U(t):=\exp\big(-itH_{0}\big{missing})$ is bounded for all time:

\|U(t)\|\leq 1,\forall t\in\mathbb{R}.

Here we used the spectral norm. This can be seen from the fact that $H_{0}$ is block diagonal, $H_{s}$ is Hermitian, and $d_{k}$ has non-negative imaginary parts.

By separating the $\mathcal{O}(\lambda)$ term in Eq. 46, we can write Eq. 44 in a perturbative form,

\partial_{t}\Gamma=-i(H_{0}\Gamma-\Gamma H_{0}^{\dagger})+\sum_{k=1}^{K}V_{k}\Gamma V_{k}^{\dagger}-i\lambda(H_{1}\Gamma-\Gamma H_{1}^{\dagger}).

(49)

In particular, we let $\Gamma^{(0)}$ be the solution of the “unperturbed” equation,

\partial_{t}\Gamma=-i(H_{0}\Gamma-\Gamma H_{0}^{\dagger})+\sum_{k=1}^{K}V_{k}\Gamma V_{k}^{\dagger}.

(50)

To further simplify notations, let $\mathcal{L}_{0}$ be the corresponding operator on the right hand side of Eq. 50. Then one can express the solution concisely as,

\Gamma^{(0)}(t)=\exp(t\mathcal{L}_{0})\Gamma^{(0)}(0).

(51)

The following Lemma provides a bound for the perturbation term in Eq. 49.

Lemma 8.

Let $\Lambda=\max_{1\leq k\leq K}\|T_{k}\|_{1}$ , and let

\Xi=H_{1}\Gamma-\Gamma H_{1}^{\dagger},

(52)

be the commutator for any Hermitian matrix $\Gamma$ . Then the trace of the first diagonal block of $\Xi$ is bounded by,

\absolutevalue{\tr\big(\Xi_{0,0}\big{missing})}\leq 4\Lambda\sum_{k>0}\|\Gamma_{0,k}\|.

(53)

For the remaining diagonal blocks, it holds that,

\sum_{k>0}\absolutevalue{\tr\big(\Xi_{k,k}\big{missing})}\leq 4\Lambda\sum_{0<j<k\leq 4K}\left\|\Gamma_{j,k}\right\|.

(54)

These estimates can be obtained by direct calculations. For instance, the first diagonal is given by,

\Xi_{0,0}=\sum_{k=1}^{K}\Big{(}\Gamma_{0,4k-3}T_{k}-T_{k}^{\dagger}\Gamma_{4k-3,0}+\Gamma_{0,4k-2}T_{k}^{\dagger}-T_{k}\Gamma_{4k-2,0}\Big{)}.

(55)

Thus the bound Eq. 53 follows from the triangle inequalities, together with the von Neumann’s trace inequality. The important observation is that the trace of the perturbation term in Eq. 49 is only controlled by the norms of the off-diagonal blocks of $\Gamma.$

We now show that the “unperturbed part” in Eq. 49, i.e., the solution of the GQME in Eq. 44 when $\lambda=0,$ has bounded solutions. The $\lambda>0$ case can then be handled using a perturbation technique.

Lemma 9.

Assume that the imaginary parts of $d_{k}$ are non-negative, i.e., $\nu_{k}\geq 0$ for all $k$ . Assume $\lambda=0$ . The solution of the GQME in Eq. 44 is denoted by $\Gamma^{(0)}(t)$ (also see Eq. 50). Then the following statements hold.

1.

The first diagonal block is given by,

$\rho(t)=U_{S}(t)\rho_{S}(0)U_{S}(t)^{\dagger},\quad U_{S}(t):=\exp(-itH_{s}).$
2.

If $\Gamma^{(0)}(0)$ is block diagonal, then $\Gamma^{(0)}(t)$ remains block diagonal.
3.

If $\Gamma^{(0)}(0)$ is given by Eq. 45, then, the trace of $\Gamma^{(0)}(t)$ is

$\displaystyle\tr\big(\Gamma^{(0)}(t)\big{missing})=2K+1+2\sum_{k}e^{-4\nu_{k}t}.$ (56)
4.

The solution $\Gamma^{(0)}(t)$ of Eq. 50 is bounded for general initial conditions. Namely, the exists a constant $c$ , independent of $t$ , such that,

$\|\Gamma^{(0)}(t)\|_{\infty}\leq c\|\Gamma^{(0)}(0)\|_{\infty},\quad\forall t\in\mathbb{R}_{+}.$

We included the proof in Appendix B.

The GQME in Eq. 44 is considered as a route to obtain an approximation to the density matrix $\rho_{S}(t)$ from Eq. 27. Assuming that the representation of the bath correlation function in Eq. 5 is exact, we can show that the error associated with the approximation of $\rho_{S}(t)$ by the GQME in Eq. 44 is ${\mathcal{O}\left(\lambda^{3}\right)}$ .

Theorem 10.

Let $\rho_{S}(t)$ be the density matrix from the full quantum model in Eq. 27 with bath correlation given by Eq. 5. In addition, let $\widehat{\rho}_{S}(t)$ be the first diagonal block of the density matrix $\Gamma(t)$ from the GQME in Eq. 44. Then,

\rho_{S}(t)=\widehat{\rho}_{S}(t)+{\mathcal{O}\left(\lambda^{3}\right)}.

(57)

This asymptotic error is consistent with that of the NMSSE in Eq. 30.

In the computation, we work with the GQME in Eq. 44, and the next theorem provides some bounds on the solution $\Gamma(t)$ .

Theorem 11.

For any $t>0$ , the density matrix $\Gamma(t)$ from the GQME in Eq. 44 is positive semidefinite: $\Gamma(t)\geq 0$ , and it has the following properties.

(i)

The norm of the density matrix follows the bound,

$\|\Gamma(t)\|_{\infty}\leq\|\Gamma(0)\|_{\infty}\exp(2\lambda C\|H_{1}\|t).$ (58)

The constant $C$ is the same as that in the Lemma 4.
(ii)

The norms of the off-diagonal blocks of $\Gamma(t)$ is of order $\lambda$ .

(iii)

For any initial condition $\Gamma(0),$ not necessarily positive, the trace of $\Gamma(t)$ is bounded as,

\mid\tr\big(\Gamma(t)\big{missing})\mid\leq 3K\mid\tr\big(\Gamma_{1,1}(0)\big{missing})\mid+\sum_{k>1}\mid\tr\big(\Gamma_{k,k}(0)\big{missing})\mid.

(59)

Corollary 12.

Fix $T>0$ . The trace of $\Gamma(t)$ from Eq. 44 for $t\in[0,T]$ is bounded as follows,

\tr\left(\Gamma(t)\right)=2K+1+2\sum_{k=1}^{K}e^{-4\nu_{k}t}+\mathcal{O}(\lambda^{2}).

(60)

More importantly, the trace of $\widetilde{\rho}_{S}(t)$ , as the first block diagonal of $\Gamma(t)$ , is bounded by,

\tr\left(\widetilde{\rho}_{S}(t)\right)=\tr\left(\Gamma_{0,0}(t)\right)=1+\mathcal{O}(\lambda^{3}).

(61)

4 Quantum algorithm for simulating generalized quantum master equations

Before presenting the quantum algorithm for this problem, we first prove some useful technical results that will be used in the analysis of our quantum algorithm.

Lemma 13 (Initial state preparation).

Given a copy of any initial state $\ket{\psi}$ of the system, a normalized version of $\Gamma_{0}$ , as defined in Eq. 45, can be prepared using $\mathcal{O}(K)$ 1- and 2-qubit gates.

Proof.

Let $K^{\prime}\leq 2K+1$ be the smallest integer larger than $K$ such that $(4K^{\prime}+1)$ is some power of 2. We initialize a $\log(4K^{\prime}+1)$ -qubit state which is a normalized version of

\displaystyle\ket{0}+\sum_{k=1}^{K}\Big{(}\ket{4k-2}+\ket{4k}\Big{)}.

(62)

This can be implemented using $\mathcal{O}(K)$ 1- and 2-qubit gates. Append an additional register that is initialized to $\ket{0}$ and use $\mathrm{CNOT}$ gates to “copy” the basis states to the second register. We have the normalized version of

\displaystyle\ket{0}\!\ket{0}+\sum_{k=1}^{K}\Big{(}\ket{4k-2}\!\ket{4k-2}+\ket{4k}\!\ket{4k}\Big{)}.

(63)

Then, simply by appending the original state $\rho_{S}$ of the system and tracing out the second register, we obtain the normalized version of $\Gamma_{0}$ . ∎

Constructing block-encodings. In the next two lemmas, we show how to obtain block-encodings of $H$ and $V_{j}$ in Eq. 44 for $j\in[K]$ using block-encodings of $H_{S}$ and $S_{\beta}$ for $\beta\in[M]$ .

Lemma 14.

Let $H$ be defined in Eq. 41, and recall $\norm{\theta}_{1}=\sum_{k=1}^{K}\theta_{k}$ . Given the access to an $(\alpha,a,\epsilon)$ -block-encoding $U_{H_{S}}$ of $H_{S}$ , an $(\alpha,a,\epsilon)$ -block-encoding $U_{S_{\beta}}$ for each $\beta\in[M]$ , an $(\alpha^{\prime},a^{\prime},\epsilon^{\prime})$ -block-encoding of $I-i\delta H$ , where

$\displaystyle\alpha^{\prime}$	$\displaystyle\leq 1+\delta\left(\alpha+d_{\max}+4\lambda\alpha\sqrt{M}\norm{\theta}_{1}\right),$	(64)
$\displaystyle a^{\prime}$	$\displaystyle=a+\mathcal{O}(\log(MK)),$	(65)
$\displaystyle\epsilon^{\prime}$	$\displaystyle=\alpha^{\prime}\epsilon.$	(66)

can be implemented with one invocation of $U_{H_{S}}$ and $2K$ invocations to each of $U_{S_{\beta}}$ together with $\mathcal{O}(MK)$ 1- and 2-qubit gates.

Proof.

Combining Eqs. 41 and 36, we can write $H$ as

\displaystyle\begin{aligned} H&=I_{A}\otimes H_{S}+H_{A}\otimes I_{S}+i\lambda\sum_{k=1}^{K}\sum_{\beta=1}^{M}\theta_{k}\innerproduct{Q_{k}}{\beta}D_{k}\otimes S_{\beta}\\ &\quad\quad+i\lambda\sum_{k=1}^{K}\sum_{\beta=1}^{M}\theta_{k}\innerproduct{\beta}{Q_{k}}E_{k}\otimes S_{\beta}^{{\dagger}}.\end{aligned}

(67)

We first show that it is easy to construct a block-encoding of the tensor product of two block-encodings. Let $U_{A}$ be an $(\alpha_{1},a,\epsilon)$ -block-encoding of $A$ and $U_{B}$ be an $(\alpha_{2},b,\epsilon)$ -encoding of $B$ , then it is straightforward to see that the unitary $U_{A}\otimes U_{B}$ together with $\mathcal{O}(a+b)$ swap gates is an $(\alpha_{1}\alpha_{2},a+b,\epsilon)$ -block-encoding of $A\otimes B$ . As a result, an $(\alpha,a,\epsilon)$ -block-encoding of $I_{A}\otimes H_{S}$ can be implemented. Since $H_{A}$ is diagonal, it is easy to implement a $(d_{\max},\mathcal{O}(\log K),0)$ -block-encoding of $H_{A}$ . Hence, a $(d_{\max},\mathcal{O}(\log K),0)$ -block-encoding of $H_{A}\otimes I_{S}$ can be implemented.

Since $D_{k}$ is $2$ -sparse, we use Lemma 4 to implement an $(2,\mathcal{O}(\log(K)),0)$ -block-encoding $U_{D_{k}}$ of $D_{k}$ . Note that this is an exact block-encoding because each nonzero entry of $D_{k}$ is 1. We also need to implement the sparse-access oracles for $D_{k}$ . Since it only has constant nonzero entries, this can be done by a small classical circuit and turning it into a quantum circuit with $\mathcal{O}(\log K)$ 1- and 2-qubit gates. Using the observation on tensor products of block-encodings, we can implement a $(2\alpha,a+\mathcal{O}(\log K),\epsilon)$ -block-encoding of $D_{k}\otimes S_{\beta}$ given an $(\alpha,a,\epsilon)$ -block-encoding $U_{S_{\beta}}$ of $S_{\beta}$ . $E_{k}\otimes S_{\beta}^{{\dagger}}$ can be implemented similarly.

Note that $I-i\delta H$ contains $3+2KM$ terms. They are $I_{A}\otimes I_{S}$ , $-\delta I_{A}\otimes H_{S}$ , $-\delta H_{A}\otimes I_{S}$ , $-i\delta\lambda\theta_{k}\innerproduct{Q_{k}}{\beta}D_{k}\otimes S_{\beta}$ , and $-i\delta\lambda\theta_{k}\innerproduct{Q_{k}}{\beta}E_{k}\otimes S_{\beta}^{{\dagger}}$ for each $k,\beta$ . The linear combination of normalizing factors is

		$\displaystyle 1+\delta\left(\alpha+d_{\max}+2\lambda\sum_{k=1}^{K}\sum_{\beta=1}^{M}\theta_{k}\lvert\innerproduct{Q_{k}}{\beta}\rvert 2\alpha\right)$		(68)
		$\displaystyle\quad\quad\leq 1+\delta\left(\alpha+d_{\max}+4\lambda\alpha\sqrt{M}\norm{\theta}_{1}\right),$		(69)

where we have used the Cauchy–Schwarz inequality. As a result, we can use Lemma 5 to construct a block-encoding of $I-i\delta H$ with the desired parameters. This construction uses $2K$ invocations to each of $U_{S_{\beta}}$ , $2M$ invocations to the block-encodings of each $D_{k}$ and $E_{k}$ (using $\mathcal{O}(M\log K)$ total 1- and 2-qubit gates), and one invocation to $U_{H_{S}}$ . The additional 1- and 2-qubit gates are used for the state preparation in Lemma 5 and this cost is $\mathcal{O}(MK)$ , which dominates the gate cost.

∎

Lemma 15.

Given $d_{k}$ ’s as in Eq. 5, a $(d_{\max},\mathcal{O}(\log(K)),\epsilon)$ -block-encoding for $V_{k}$ for $k\in[K]$ can be constructed using $\mathcal{O}(\log K+\mathrm{polylog}(1/\epsilon))$ 1- and 2-qubit gates.

Proof.

Each $V_{k}$ is 1-sparse, and it is straightforward to implement the sparse-access oracles specified in Eq. 13 and Eq. 14. Then we use Lemma 4 to implement a $(1,\mathcal{O}(\log(K)),\epsilon)$ -block-encoding for $V_{k}/d_{\max}$ , which implies a $(d_{\max},\mathcal{O}(\log(K)),\epsilon)$ -block-encoding for $V_{k}$ . ∎

Infinitesimal approximation by completely-positive maps. An important step of our quantum algorithm is to use the following superoperator to approximate $e^{\mathcal{K}\delta}$ when $\delta$ is small.

\displaystyle\mathcal{M}_{\delta}(X)=A_{0}XA_{0}^{{\dagger}}+\sum_{k=1}^{K}A_{k}XA_{k}^{{\dagger}},

(70)

where

\displaystyle A_{0}=I-i\delta H,\text{ and }\quad A_{k}=\sqrt{\delta}V_{k}\text{ for all $k\in[K]$}.

(71)

We use the following lemma, which is proved in Appendix F, to bound the error of this approximation:

Lemma 16.

Let $\mathcal{M}_{\delta}$ be a superoperator defined in Eq. 70, and let $\mathcal{K}$ be as defined in Eq. 44. Define $\Lambda:=\norm{H}+\sum_{k=1}^{K}\norm{V_{k}}^{2}$ . Then it holds that

\displaystyle\norm{\mathcal{M}_{\delta}-e^{\mathcal{K}\delta}}_{\diamond}\leq 5(\delta\Lambda)^{2}.

(72)

Oblivious amplitude amplification for isometries. In our algorithm, we need to apply amplitude amplification to boost the success probability. However, in our context, the underlying operator is an isometry instead of a unitary (i.e., part of the input is restricted to be some special state). We use oblivious amplitude amplification for isometries, which was first introduced in [CW17] to achieve this. In [CW17], only a special case where the initial success probability is exactly $1/4$ was considered. Here, we give a more general version.

Lemma 17.

For any $a,b\in\mathbb{N}_{+}$ , let $\ket*{\widehat{0}}:=\ket{0}^{\otimes a}$ and $\ket*{\widehat{\mu}}:=\ket{\mu}^{\otimes b}$ for an arbitrary state $\ket{\mu}$ . For any $n$ -qubit state $\ket{\psi}$ , define $\ket*{\widehat{\psi}}:=\ket*{\widehat{0}}\ket*{\widehat{\mu}}\ket{\psi}$ . Let the target state $\ket*{\widehat{\phi}}$ be defined as

\displaystyle\ket*{\widehat{\phi}}:=\ket*{\widehat{0}}\ket{\phi},

(73)

where $\ket{\phi}$ is a $(b+n)$ -qubit state. Let $P_{0}:=\outerproduct*{\widehat{0}}{\widehat{0}}\otimes I_{2^{b}}\otimes I_{2^{n}}$ and $P_{1}:=\outerproduct*{\widehat{0}}{\widehat{0}}\otimes\outerproduct{\widehat{\mu}}{\widehat{\mu}}\otimes I_{2^{n}}$ be two projectors. Suppose there exists an operator $W$ such that

\displaystyle W\ket*{\widehat{\psi}}=\sin\theta\ket*{\widehat{\phi}}+\cos\theta\ket*{\widehat{\phi}^{\bot}},

(74)

for some $\theta\in[0,\pi/2]$ with $\ket*{\widehat{\phi}^{\bot}}$ satisfying $P_{0}\ket*{\widehat{\phi}^{\bot}}=0$ . Then it holds that

\displaystyle-\!W\!\left(I\!-\!2P_{1}\right)W^{{\dagger}}\left(I\!-\!2P_{0}\right)(\sin\gamma\ket*{\widehat{\phi}}\!+\!\cos\gamma\ket*{\widehat{\phi}^{\bot}})\!=\!\sin(\gamma\!+\!2\theta)\ket*{\widehat{\phi}}\!+\!\cos(\gamma\!+\!2\theta)\ket*{\widehat{\phi}^{\bot}},

(75)

for all $\gamma\in[0,\pi/2]$ .

Proof.

Let $\ket*{\widehat{\psi}^{\bot}}$ be a state satisfying

\displaystyle W\ket*{\widehat{\psi}^{\bot}}=\cos\theta\ket*{\widehat{\phi}}-\sin\theta\ket*{\widehat{\phi}^{\bot}}.

(76)

It is useful to have the following facts

	$\displaystyle W^{{\dagger}}\ket*{\widehat{\phi}}$	$\displaystyle=\sin\theta\ket{\widehat{\psi}}+\cos\theta\ket{\widehat{\psi}^{{\dagger}}}\quad\text{ and}$		(77)
	$\displaystyle W^{{\dagger}}\ket*{\widehat{\phi}^{\bot}}$	$\displaystyle=\cos\theta\ket{\widehat{\psi}}-\sin\theta\ket{\widehat{\psi}^{{\dagger}}},$		(78)

which can be obtained from Eq. 74 and Eq. 76.

We first show that $P_{1}\ket*{\widehat{\psi}^{\bot}}=0$ . To see this, we define an operator

\displaystyle Q=(\bra*{\widehat{0}}\bra*{\widehat{\mu}}\otimes I)W^{{\dagger}}P_{0}W(\ket*{\widehat{0}}\ket*{\widehat{\mu}}\otimes I).

(79)

For any $\ket{\psi}$ , we have

	$\displaystyle\bra{\psi}Q\ket{\psi}$	$\displaystyle=\norm{P_{0}W\left(\ket{\widehat{0}}\ket{\widehat{\mu}}\ket{\psi}\right)}^{2}$		(80)
		$\displaystyle=\norm{P_{0}\left(\sin\theta\ket{\widehat{\phi}}+\cos\theta\ket{\widehat{\phi}^{\bot}}\right)}^{2}=\norm{\sin\theta\ket*{\widehat{\phi}}}^{2}=\sin^{2}\theta.$		(80)

Hence, all the eigenvalues of $Q$ are $\sin^{2}\theta$ , so we can write

\displaystyle Q=\sin^{2}\theta I.

(81)

Now, consider any state $\ket{\psi}$ :

$\displaystyle Q\ket{\psi}$	$\displaystyle=(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)W^{{\dagger}}P_{0}W(\ket{\widehat{0}}\ket{\widehat{\mu}}\otimes\ket{\psi})=\sin\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)W^{{\dagger}}\ket*{\widehat{\phi}}$	(82)
	$\displaystyle=\sin\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)\left(\sin\theta\ket{\widehat{\psi}}+\cos\theta\ket{\widehat{\psi}^{\bot}}\right)$
	$\displaystyle=\sin^{2}\theta\ket{\psi}+\sin\theta\cos\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)\ket*{\widehat{\psi}^{\bot}},$

where the third equality follows from Eq. 74 and Eq. 76. On the other hand, by Eq. 81, we have

\displaystyle Q\ket{\psi}=\sin^{2}\theta\ket{\psi}.

(83)

In light of Eq. 82 and Eq. 83, we must have,

\displaystyle\left(\bra*{\widehat{0}}\!\bra*{\widehat{\mu}}\otimes I\right)\ket*{\widehat{\psi}^{\bot}}=0,

(84)

which implies that $P_{1}\ket*{\widehat{\psi}^{\bot}}=0$ .

We now proceed to analyze the result of applying $W(I-2P_{1})W^{{\dagger}}(I-2P_{0})$ on $\sin\gamma\ket*{\widehat{\phi}}+\cos\gamma\ket*{\widehat{\phi}^{\bot}}$ ,

	$\displaystyle W(I-2P_{1})W^{{\dagger}}(I-2P_{0})\left(\sin\gamma\ket{\widehat{\phi}}+\cos\gamma\ket{\widehat{\phi}^{\bot}}\right)$		(85)
	$\displaystyle=\left(I-2P_{0}-2WP_{1}W^{{\dagger}}+4WP_{1}W^{{\dagger}}P_{0}\right)\left(\sin\gamma\ket{\widehat{\phi}}+\cos\gamma\ket{\widehat{\phi}^{\bot}}\right)$		(86)
	$\displaystyle=\sin(\gamma+2\theta)\ket{\widehat{\phi}}+\cos(\gamma+2\theta)\ket{\widehat{\phi}^{\bot}}.$		(87)

∎

Proof of the main theorem. Now, we have all the tools for proving the quantum algorithm for simulating 1.

Theorem 18.

Suppose we are given a block-encoding $U_{H_{S}}$ of $H_{S}$ , a block-encodings $U_{S_{\alpha}}$ of $S_{\alpha}$ (for $\alpha\in[M]$ ), $\theta_{k}$ , $\ket{Q_{k}}$ (all entries), $\lambda$ , and $d_{k}$ for $k\in[K]$ as in Eq. 5. There exists a quantum algorithm that solves 1. Let $\tau=t(\alpha(1+\lambda\sqrt{M}\norm{\theta}_{1})+Kd_{\max}^{2})$ , where $\norm{\theta}_{1}=\sum_{k=1}^{K}\theta_{k}$ . This quantum algorithm uses

\displaystyle\mathcal{O}\left(\tau\sqrt{K}\frac{\log(\tau/\epsilon)}{\log\log(\tau/\epsilon)}\right),

(88)

queries to $U_{H_{S}}$ and $U_{S_{a}}$ , and

\displaystyle\mathcal{O}\left(\frac{\tau\sqrt{K}\log^{2}(\tau/\epsilon)}{\log\log(\tau/\epsilon)}+\tau K^{2.5}\right),

(89)

additional 1- and 2-qubit gates.

We first use Lemma 13 to prepare a normalized version of $\Gamma_{0}$ as in Eq. 45. To simulate the dynamics in Eq. 44, first note that the solution to Eq. 44 is $e^{\mathcal{K}t}$ . We use the superoperator $\mathcal{M}_{\delta}$ defined in Eq. 70 and Eq. 71 to approximate $e^{\delta\mathcal{K}}$ for a small step $\delta$ . By Lemma 16, we know that the approximation error is at most $5(\delta\Lambda)^{2}$ , where $\Lambda=\norm{H}+\sum_{k=1}^{K}\norm{V_{k}}^{2}$ .

To use Lemma 6 to implement $\mathcal{M}_{\delta}$ , we first need to implement the block-encoding of $A_{j}$ for $j\in\{0,\ldots,K-1\}$ . Using Lemma 14, a $(\alpha^{\prime},a^{\prime},\epsilon^{\prime})$ -block-encoding of $A_{0}$ can be implemented with $\alpha^{\prime}\leq 1+\delta(\alpha+d_{\max}+4\lambda\alpha\sqrt{M}\norm{\theta}_{1}))$ , $a^{\prime}=a+\mathcal{O}(\log(MK))$ , and $\epsilon^{\prime}=\alpha^{\prime}\epsilon$ . Also, a $(\sqrt{\delta}d_{\max},\mathcal{O}(\log(K)),\epsilon)$ -block-encoding of $A_{j}$ can be implemented for each $j\in[K]$ using $\mathcal{O}(\log K+\mathrm{polylog}(1/\epsilon))$ 1- and 2-qubit gates by Lemma 4. Now we use Lemma 6 to obtain the unnormalized state $\sum_{j=1}^{m}\ket{j}A_{j}\ket{\psi}$ with the “success probability” parameter

	$\displaystyle\geq\frac{1}{\left(1+\delta(\alpha+d_{\max}+4\lambda\alpha\sqrt{M}\norm{\theta}_{1})\right)^{2}+\sum_{j=1}^{K}\delta d_{\max}^{2}}$		(90)
	$\displaystyle=\frac{1}{(1+\delta(2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2})+\delta^{2}(\alpha+d_{\max}+4\alpha\lambda\sqrt{M}\norm{\theta}_{1})^{2}}$		(91)
	$\displaystyle=1-\delta(2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2})+\mathcal{O}(\delta^{2}(\alpha+d_{\max}+4\alpha\lambda\sqrt{M}\norm{\theta}_{1})^{2}).$		(92)

Setting the step size,

\displaystyle\delta=\Theta\left(\frac{1}{r(2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2})}\right),

(93)

and repeating the above procedure $r$ times, this success probability parameter becomes $\mathcal{O}(1)$ .

Now for the approximation error, we observe that

		$\displaystyle\norm{\mathcal{M}_{\delta}-e^{\delta\mathcal{K}}}_{\diamond}$		(94)
		$\displaystyle\leq 5\delta^{2}\Lambda^{2}=\mathcal{O}\left(\frac{\Lambda^{2}}{r^{2}(2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2})^{2}}\right)\leq\mathcal{O}\left(\frac{1}{r^{2}}\right).$		(94)

By applying the approximation $r$ times, we have

\displaystyle\norm{\mathcal{M}_{\delta}^{r}-e^{t\mathcal{K}}}_{\diamond}\leq\mathcal{O}\left(\frac{1}{r}\right),

(95)

for evolution $t=r\delta=\Theta\left(\frac{1}{2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M\norm{\theta}_{1}}+Kd_{\max}^{2}}\right)$ . The parameter $r$ can be chosen larger enough so that this error is at most $\epsilon$ .

Further, conditioned on that we have obtained an approximation $\widetilde{\Gamma}_{t}$ of $\Gamma_{t}$ , we can extract an approximation of $\outerproduct{\psi}{\psi}$ by measuring the first $\log(4K+1)$ qubits of $\widetilde{\Gamma}_{t}$ and post-select the outcome being 0. According to Corollary 12, this probability of the outcome being 0 is $\Omega(1/(2K+1))$ . Now, the success probability is $\Omega(1/K)$ . It follows from Lemma 17 that using $\mathcal{O}(\sqrt{K})$ iterations of oblivious amplitude amplification for isometries, we obtain the desired state.

Now, the total evolution time we have simulated so far is $r\delta=\Theta\left(\frac{1}{2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2}}\right)$ , and the total cost in terms of queries to $U_{H_{S}}$ and $U_{S_{\beta}}$ is $\mathcal{O}(r)$ . In the following, we show how to achieve poly-logarithmic cost. Note that when applying Lemma 5 to construct a block-encoding of $A_{0}=I-i\delta H$ and applying Lemma 6 to implement the superoperator specified by Kraus operators $A_{0},A_{1},\ldots,A_{K}$ , the first register (containing $\mathcal{O}(\log K)$ qubits) is used for the $\ket{\mu}$ state in Lemma 6, and the second register (containing $\mathcal{O}(\log K)$ qubits) is used for state $B\ket{0}$ in Lemma 5. The coefficients of the state in the two registers are concentrated to $\ket{0}\!\ket{0}$ , which corresponds to $I$ (nothing to implement). More specifically, recall that the parameter $s_{0}$ in Lemma 6 is the block-encoding normalization factor for $I-i\delta H$ , which can be written as $s_{0}=\sum_{j=1}^{2KM+3}y_{i}\alpha_{j}$ , where $y_{j}$ ’s and $\alpha_{j}$ ’s are the parameters as used in Lemma 5. For convenience, let $\widehat{\alpha}$ be upper bounded by $\widehat{\alpha}\leq\alpha+d_{\max}+4\alpha\lambda\sqrt{M}\norm{\theta}_{1}$ , so we can write $\sum_{j=1}^{2KM+3}y_{i}\alpha_{j}=1+\delta\widehat{\alpha}$ according to Lemma 14. As a result, the amplitude for $\ket{0}\ket{0}$ is

	$\displaystyle\frac{s_{0}}{\sqrt{\sum_{j=0}^{K}s_{j}^{2}}}\cdot\frac{1}{\sqrt{\sum_{j=1}^{2KM+3}y_{j}\alpha_{j}}}=\sqrt{\frac{\sum_{j=1}^{2KM+3}y_{j}\alpha_{j}}{\sum_{j=0}^{K}s_{j}^{2}}}$		(96)
	$\displaystyle=\sqrt{\frac{1+\delta\widehat{\alpha}}{(1+\delta\widehat{\alpha})^{2}+Kd_{\max}^{2}\delta}}$		(97)
	$\displaystyle=\sqrt{\frac{1+\delta(\widehat{\alpha})}{1+2\delta\widehat{\alpha}+\delta Kd_{\max}^{2}\!+\!\delta^{2}\widehat{\alpha}^{2}}}$		(98)
	$\displaystyle=\sqrt{1-\delta(\widehat{\alpha}+Kd_{\max}^{2})+\Theta(\delta^{2}(2\widehat{\alpha}+Kd_{\max}^{2})^{2})}.$		(99)

Therefore, the probability that the first two registers are not measured $0,0$ is proportional to³³3Note that we use the term “proportional to” because the actual probability is normalized according to the trace ratio of the first block and the whole matrix of $\Gamma$ , which incurs a factor $O(1/K)$ .

\displaystyle\delta\!\left(\widehat{\alpha}+Kd_{\max}^{2}\right)+\Theta\left(\delta^{2}(2\widehat{\alpha}+Kd_{\max}^{2})^{2}\right)=\mathcal{O}\left(\frac{1}{r}\right).

(100)

We apply the techniques used in [CW17] (first introduced in [BCC⁺14]) to bypass the state preparation for $\ket{\mu}$ in Lemma 5 and for $B\ket{0}$ in Lemma 5. Instead, we use a state with Hamming weight at most

\displaystyle h=\mathcal{O}\left(\frac{\log(1/\epsilon)}{\log\log(1/\epsilon)}\right),

(101)

to approximate the state after the state preparation procedure while causing error at most $\epsilon$ . Following similar analysis as in [CW17], the number of 1- and 2-qubit gates for implementing this compressed encoding procedure is $\mathcal{O}(\log(1/\epsilon)h+K^{2})$ . Now, the number of queries to $U_{H_{S}}$ and $U_{S_{\beta}}$ becomes

\displaystyle\mathcal{O}\left(\frac{\sqrt{K}\log(1/\epsilon)}{\log\log(1/\epsilon)}\right).

(102)

For arbitrary simulation time $t$ , we repeat this $\mathcal{O}(t(2\alpha+2d_{\max}+8\alpha\lambda\sqrt{M}\norm{\theta}_{1}+Kd_{\max}^{2}))=\mathcal{O}(t(\alpha(1+\lambda\sqrt{M}\norm{\theta}_{1})+Kd_{\max}^{2}))$ times where each segment has the error parameter $\epsilon/(t(\alpha(1+\lambda\sqrt{M}\norm{\theta}_{1})+Kd_{\max}^{2}))$ , which has total cost

\displaystyle\mathcal{O}\left(\tau\sqrt{K}\frac{\log(\tau/\epsilon)}{\log\log(\tau/\epsilon)}\right).

(103)

Here we have $\tau=t(\alpha(1+\lambda\sqrt{M}\norm{\theta}_{1})+Kd_{\max}^{2})$ .

The additional 1- and 2-qubit gates are used in the compressed encoding procedure and the reflections in the oblivious amplitude amplification for isometries. This is bounded by

\displaystyle\mathcal{O}\left(\tau\sqrt{K}\left(\log(\tau/\epsilon)h+K^{2}\right)\right)=\mathcal{O}\left(\frac{\tau\sqrt{K}\log^{2}(\tau/\epsilon)}{\log\log(\tau/\epsilon)}+\tau K^{2.5}\right).

(104)

5 Acknowledgments

CW thanks Yudong Cao and Peter D. Johnson for helpful discussions on the HEOM approach for modeling non-Markovian open quantum systems. XL’s research is supported by the National Science Foundation Grants DMS-2111221.

References

[AM21] Alexey Alliluev and Denis Makarov. Dynamics of a nonlinear quantum oscillator under non-Markovian pumping. arXiv preprint arXiv:2110.07914, 2021.
[BACS07] Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders. Efficient quantum algorithms for simulating sparse Hamiltonians. Communications in Mathematical Physics, 270(2):359–371, 2007.
[BCC⁺14] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma. Exponential improvement in precision for simulating sparse Hamiltonians. In Proceedings of the 46th annual ACM symposium on Theory of computing (STOC 2014), pages 283–292, 2014.
[BCC⁺15] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma. Simulating Hamiltonian dynamics with a truncated Taylor series. Physical Review Letters, 114(9), 2015.
[BCC⁺17] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma. Exponential improvement in precision for simulating sparse Hamiltonians. Forum of Mathematics, Sigma, 5, 2017.
[BCG14] Dominic W Berry, Richard Cleve, and Sevag Gharibian. Gate-efficient discrete simulations of continuous-time quantum query algorithms. Quantum Information & Computation, 14(1-2):1–30, 2014.
[BCK15] Dominic W Berry, Andrew M Childs, and Robin Kothari. Hamiltonian simulation with nearly optimal dependence on all parameters. In Proceedings of the 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS 2015), 2015.
[BD12] Robert Biele and Roberto D’Agosta. A stochastic approach to open quantum systems. Journal of Physics: Condensed Matter, 24(27):273201, 2012.
[BLP09] Heinz-Peter Breuer, Elsi-Mari Laine, and Jyrki Piilo. Measure for the degree of non-Markovian behavior of quantum processes in open systems. Physical Review Letters, 103(21):210401, 2009.
[BLPV16] Heinz-Peter Breuer, Elsi-Mari Laine, Jyrki Piilo, and Bassano Vacchini. Non-Markovian dynamics in open quantum systems. Reviews of Modern Physics, 88(2), 2016.
[BP02] Heinz-Peter Breuer and Francesco Petruccione. The theory of open quantum systems. Oxford University Press on Demand, Oxford, 2002.
[CGJ19] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery. The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), 2019.
[CK10] Dariusz Chruściśki and Andrzej Kossakowski. Non-Markovian Quantum Dynamics: Local versus Nonlocal. Physical Review Letters, 104(7):070406, February 2010.
[CL17a] Yu Cao and Jianfeng Lu. Lindblad equation and its semiclassical limit of the Anderson-Holstein model. Journal of Mathematical Physics, 58(12):122105, 2017.
[CL17b] Andrew M Childs and Tongyang Li. Efficient simulation of sparse Markovian quantum dynamics. Quantum Information & Computation, 17(11&12):0901–0947, 2017.
[CW17] Richard Cleve and Chunhao Wang. Efficient quantum algorithms for simulating Lindblad evolution. In 44th International Colloquium on Automata, Languages, and Programming, (ICALP 2017), pages 17:1–17:14, 2017.
[DS97] Lajos Diósi and Walter T Strunz. The non-Markovian stochastic Schrödinger equation for open systems. Physics Letters A, 235(6):569–573, 1997.
[dVA17] Inés de Vega and Daniel Alonso. Dynamics of non-Markovian open quantum systems. Reviews of Modern Physics, 89(1), 2017.
[GN99] Pierre Gaspard and Masataka Nagaoka. Non-Markovian stochastic Schrödinger equation. The Journal of Chemical Physics, 111(13):5676–5690, 1999.
[GN05] Michael Galperin and Abraham Nitzan. Current-induced light emission and light-induced current in molecular-tunneling junctions. Physical Review Letters, 95(20):206802, 2005.
[GSLW19] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC 2019), pages 193–204, 2019.
[GTP⁺15] S Gröblacher, A Trubarov, N Prigge, GD Cole, M Aspelmeyer, and J Eisert. Observation of non-Markovian micromechanical Brownian motion. Nature Communications, 6:7606, 2015.
[IT05] Akihito Ishizaki and Yoshitaka Tanimura. Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach. Journal of the Physical Society of Japan, 74(12):3131–3134, 2005.
[JZY08] Jinshuang Jin, Xiao Zheng, and YiJing Yan. Exact dynamics of dissipative electronic systems and quantum transport: Hierarchical equations of motion approach. The Journal of Chemical Physics, 128(23):234703, 2008.
[KBG⁺11] Martin Kliesch, Thomas Barthel, Christian Gogolin, Michael J Kastoryano, and Jens Eisert. Dissipative quantum Church-Turing theorem. Physical Review Letters, 107(12), 2011.
[KMCWM16] Aaron Kelly, Andrés Montoya-Castillo, Lu Wang, and Thomas E Markland. Generalized quantum master equations in and out of equilibrium: When can one win? The Journal of Chemical Physics, 144(18):184105, 2016.
[KP13] Peter E Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations, volume 23. Springer Science & Business Media, Berlin, 2013.
[Kra71] Karl Kraus. General state changes in quantum theory. Annals of Physics, 64(2):311–335, 1971.
[Kub66] Rep Kubo. The fluctuation-dissipation theorem. Reports on Progress in Physics, 29(1):255, 1966.
[LBW01] Daniel A Lidar, Zsolt Bihary, and K Birgitta Whaley. From completely positive maps to the quantum Markovian semigroup master equation. Chemical Physics, 268(1-3):35–53, 2001.
[LC17] Guang Hao Low and Isaac L Chuang. Optimal Hamiltonian simulation by quantum signal processing. Physical Review Letters, 118(1), 2017.
[LC19] Guang Hao Low and Isaac L Chuang. Hamiltonian simulation by qubitization. Quantum, 3:163, 2019.
[LCW98] Daniel A Lidar, Isaac L Chuang, and K Birgitta Whaley. Decoherence-free subspaces for quantum computation. Physical Review Letters, 81(12):2594, 1998.
[LFS12] Shunlong Luo, Shuangshuang Fu, and Hongting Song. Quantifying non-Markovianity via correlations. Physical Review A, 86(4):044101, 2012.
[Li21] Xiantao Li. Markovian embedding procedures for non-Markovian stochastic schrödinger equations. Physics Letters A, 387:127036, 2021.
[Lin76] Goran Lindblad. On the generators of quantum dynamical semigroups. Communications in Mathematical Physics, 48(2):119–130, 1976.
[Llo96] Seth Lloyd. Universal quantum simulators. Science, 273(5278):1073–1078, 1996.
[LRA⁺20] Neill Lambert, Tarun Raheja, Shahnawaz Ahmed, Alexander Pitchford, and Franco Nori. BoFiN-HEOM: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics. arXiv preprint arXiv:2010.10806, 2020.
[MALH⁺11] Kristian Høeg Madsen, Serkan Ates, Toke Lund-Hansen, Andreas Löffler, Stephan Reitzenstein, Alfred Forchel, and Peter Lodahl. Observation of non-Markovian dynamics of a single quantum dot in a micropillar cavity. Physical Review Letters, 106(23):233601, 2011.
[MCR16] Andrés Montoya-Castillo and David R Reichman. Approximate but accurate quantum dynamics from the Mori formalism: I. Nonequilibrium dynamics. The Journal of Chemical Physics, 144(18):184104, 2016.
[MCR17] Andrés Montoya-Castillo and David R Reichman. Approximate but accurate quantum dynamics from the Mori formalism. II. Equilibrium time correlation functions. The Journal of Chemical Physics, 146(8):084110, 2017.
[MEL20] DV Makarov, AA Elistratov, and Yu E Lozovik. Non-Markovian effects in dynamics of exciton-polariton Bose condensates. Physics Letters A, 384(36):126942, 2020.
[MT99] Christoph Meier and David J. Tannor. Non-Markovian evolution of the density operator in the presence of strong laser fields. The Journal of Chemical Physics, 111(8):3365–3376, 1999.
[Nak58] Sadao Nakajima. On quantum theory of transport phenomena: steady diffusion. Progress of Theoretical Physics, 20(6):948–959, 1958.
[Pav14] Grigorios A Pavliotis. Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations, volume 60. Springer, New York, 2014.
[PMCKM19] William C. Pfalzgraff, Andrés Montoya-Castillo, Aaron Kelly, and Thomas E. Markland. Efficient construction of generalized master equation memory kernels for multi-state systems from nonadiabatic quantum-classical dynamics. The Journal of Chemical Physics, 150(24):244109, 2019.
[PRRF⁺18] Felix A. Pollock, César Rodr$́\text{I}$guez-Rosario, Thomas Frauenheim, Mauro Paternostro, and Kavan Modi. Non-Markovian quantum processes: Complete framework and efficient characterization. Physical Review A, 97(1):012127, 2018.
[PZA⁺19] Ricardo Puebla, Giorgio Zicari, Iñigo Arrazola, Enrique Solano, Mauro Paternostro, and Jorge Casanova. Spin-boson model as a simulator of non-Markovian multiphoton Jaynes-Cummings models. Symmetry, 11(5):695, 2019.
[RE14] Gerhard Ritschel and Alexander Eisfeld. Analytic representations of bath correlation functions for ohmic and superohmic spectral densities using simple poles. The Journal of Chemical Physics, 141(9):094101, 2014.
[Ris84] Hannes Risken. Fokker-Planck Equation. Springer, Berlin, 1984.
[SDG99] Walter T. Strunz, Lajos Diósi, and Nicolas Gisin. Open System Dynamics with Non–Markovian Quantum Trajectories. Physical Review Letters, 82(9):1801–1805, 1999.
[SES14] D. Suess, A. Eisfeld, and W. T. Strunz. Hierarchy of Stochastic Pure States for Open Quantum System Dynamics. Physical Review Letters, 113(15), 2014.
[SG03] Qiang Shi and Eitan Geva. A new approach to calculating the memory kernel of the generalized quantum master equation for an arbitrary system–bath coupling. The Journal of Chemical Physics, 119(23):12063–12076, 2003.
[SH04] K Shiokawa and BL Hu. Qubit decoherence and non-Markovian dynamics at low temperatures via an effective spin-boson model. Physical Review A, 70(6):062106, 2004.
[Sha04] Jiushu Shao. Decoupling quantum dissipation interaction via stochastic fields. The Journal of Chemical Physics, 120(11):5053–5056, 2004.
[SHMS⁺21] Anthony W Schlimgen, Kade Head-Marsden, LeeAnn M Sager, Prineha Narang, and David A Mazziotti. Quantum simulation of open quantum systems using a unitary decomposition of operators. Physical Review Letters, 127(27):270503, 2021.
[SHMSS⁺22] Anthony W Schlimgen, Kade Head-Marsden, LeeAnn M Sager-Smith, Prineha Narang, and David A Mazziotti. Quantum state preparation and non-unitary evolution with diagonal operators. arXiv preprint arXiv:2205.02826, 2022.
[SKK⁺18] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett. Efficient non-Markovian quantum dynamics using time-evolving matrix product operators. Nature Communications, 9(1):3322, 2018.
[SP16] I. Semina and F. Petruccione. The simulation of the non-Markovian behaviour of a two-level system. Physica A: Statistical Mechanics and its Applications, 450:395–402, 2016.
[SSS⁺16] Ryan Sweke, Mikel Sanz, Ilya Sinayskiy, Francesco Petruccione, and Enrique Solano. Digital quantum simulation of many-body non-Markovian dynamics. Physical Review A, 94(2):022317, 2016.
[Str96] Walter T Strunz. Linear quantum state diffusion for non–Markovian open quantum systems. Physics Letters A, 224(1-2):25–30, 1996.
[Tan06] Yoshitaka Tanimura. Stochastic Liouville, Langevin, Fokker–Planck, and master equation approaches to quantum dissipative systems. Journal of the Physical Society of Japan, 75(8):082001, 2006.
[Tan20] Yoshitaka Tanimura. Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM). The Journal of Chemical Physics, 153(2):020901, 2020.
[TER⁺09] Michael Thorwart, Jens Eckel, John H Reina, Peter Nalbach, and Stephan Weiss. Enhanced quantum entanglement in the non-Markovian dynamics of biomolecular excitons. Chemical Physics Letters, 478(4-6):234–237, 2009.
[VMP⁺11] Ruggero Vasile, Sabrina Maniscalco, Matteo GA Paris, Heinz-Peter Breuer, and Jyrki Piilo. Quantifying non-Markovianity of continuous-variable Gaussian dynamical maps. Physical Review A, 84(5):052118, 2011.
[WC13] Chen Wang and Qing-Hu Chen. Exact dynamics of quantum correlations of two qubits coupled to bosonic baths. New Journal of Physics, 15(10):103020, 2013.
[Zwa60] Robert Zwanzig. Ensemble method in the theory of irreversibility. The Journal of Chemical Physics, 33(5):1338–1341, 1960.

Appendix A The derivation of the extended stochastic dynamics Eq. 35

As a building block to approximate stationary Gaussian processes in the NMSSE (Eq. 30), we consider the Ornstein–Uhlenbeck (OU) type of processes [Ris84], expressed as the solution of the following linear SDEs,

i\dot{\zeta}_{k}=-d_{k}\zeta_{k}+\gamma_{k}\dot{w}_{j}(t),

(105)

for $k=1,2,\ldots,K$ , with $K\geq M$ . Here $\gamma_{k}\geq 0$ and $d_{k}$ is a complex number with positive imaginary part. In addition, each of these independent OU processes has an initial variance $\theta_{k}^{2}$ ,

\mathbb{E}[\zeta_{k}(0)^{\dagger}\zeta_{k}(0)]=\theta_{k}^{2}.

If we pick $\gamma_{k}$ , such that $\gamma_{k}^{2}=2\text{Im}(d_{k}),$ then $\zeta_{k}(t)$ is a stationary Gaussian process with correlation,

\mathbb{E}[\zeta_{j}^{\dagger}(t)\zeta_{k}(t^{\prime})]=\theta_{j}^{2}\delta_{j,k}\exp\big(-id_{j}^{*}(t-t^{\prime})\big{missing}).

(106)

We now construct an ansatz for approximating the bath correlation function. The case when $M=1$ has been thoroughly investigated in [RE14]. As a time correlation function, $C(t)$ can be expressed in terms of its power spectrum, denoted here by $\widehat{C}(\omega),$ as a Fourier integral,

C(t)=\int_{-\infty}^{\infty}\widehat{C}(\omega)e^{-i\omega t}d\omega.

(107)

Known as the power spectrum, $\widehat{C}(\omega)$ is Hermitian and positive semidefinite. Hence it can be diagonalized using the transformation:

\widehat{C}(\omega)=\sum_{j}\lambda_{j}(\omega)^{2}\ket{Q_{j}(\omega)}\bra{Q_{j}(\omega)}.

Therefore, a direct numerical approximation of Eq. 107 will certainly lead to an ansatz like Eq. 5.

Another alternative is to use a contour integral in the upper half plane, and the Cauchy residue theorem, reducing the Fourier integral in Eq. 107 to a summation over poles [RE14]. This approach will also lead to the ansatz like Eq. 5.

To ensure the approximation accuracy, in practice, the ansatz in Eq. 5 is often obtained by a least-squares approach. In addition, we will consider the poles $d_{k}$ within a cut-off frequency $d_{\text{max}}$ , i.e., $\absolutevalue{d_{k}}\leq d_{\text{max}}$ .

Next we show that we can find an approximation of $\eta(t)$ that exactly satisfies the relation in Eq. 31, where $C(t)$ is represented by Eq. 5. Specifically, we approximate the noise $\eta(t)$ using the OU processes $\zeta_{k}(t)$ from Eq. 105,

\eta_{\beta}(t)=\sum_{k=1}^{K}\innerproduct{Q_{k}}{\beta}\theta_{k}\zeta_{k}(t).

(108)

In light of Eq. 106, the approximation in Eq. 108 corresponds to an approximation of the time correlation function in Eq. 5.

For the case of a single interaction term, where $M=1$ and the matrix $C(t)$ is corresponding to a scalar function, this reduces to the standard approach using a sum of exponentials [RE14]. But Eq. 5 provides a more general scheme to handle multiple interaction terms.

Next we show how this can simplify the NMSSE in Eq. 30 when the bath correlation functions are expressed as Eq. 5. We first define the operators $T_{j}$ according to Eq. 36. This simplifies the NMSSE in Eq. 30 to,

i\partial_{t}\psi=\hat{H}_{S}\psi-i\lambda^{2}\sum_{k=1}^{K}\int_{0}^{t}T^{\dagger}_{k}e^{-i(\hat{H}_{S}+d_{k}^{*})\tau}T_{k}\psi(t-\tau)d\tau+\lambda\sum_{j=1}^{K}\zeta_{j}(t)T_{j}\psi(t).

(109)

Here the multiplication by $d_{k}^{*}$ represents the operator $d_{k}^{*}I_{S}$ , where $I_{S}$ is the identity matrix.

Eq. 109 still contains memory. But compared to the original NMSSE in Eq. 30, the correlation $C(t)$ has been broken down to exponential functions, which can be combined with the unitary operator $e^{-it\hat{H}_{S}}.$ In addition, the noise is expressed as the OU process $\zeta_{j}$ ’s, which can be treated using Itô calculus.

Next, we demonstrate how to embed the dynamics in Eq. 109 into an extended, but Markovian, dynamics. The simple observation that motivated the Markovian embedding is that a convolution integral in time can be represented as the solution of a differential equation:

f(t)=\int_{0}^{t}\exp(-a(t-\tau))g(\tau)d\tau\quad\Longrightarrow f^{\prime}=-af+g,\;f(0)=0.

(110)

Here $f$ will be regarded as an auxiliary variable, introduced to reduce the memory integral. Compared to computing the integral at every time step using direct quadrature formulas, it is much more efficient to solve the differential equation.

We now show that we can use the same idea to define auxiliary orbitals. Specifically, by inserting Eq. 5 in the SSE in Eq. 30, and by letting, for $k=1,2,\ldots,K,$

\chi_{k}^{\text{I}}(t)=\int_{0}^{t}e^{-i(\hat{H}_{S}+d_{k}^{*})\tau}T_{k}\psi(t-\tau)d\tau,

we arrive at the equation,

i\partial_{t}\chi_{k}^{\text{I}}=(H_{S}+d_{k}^{*})\chi_{k}^{\text{I}}+i\lambda T_{k}\psi(t).

(111)

This simplifies Eq. 109 to,

i\partial_{t}\psi=\hat{H}_{S}\psi-i\lambda^{2}\sum_{j=1}^{K}T_{j}^{\dagger}\chi^{\text{I}}_{j}+\lambda\sum_{j=1}^{K}T_{j}\psi(t)\zeta_{j}(t).

(112)

To incorporate the noise term, we define,

\chi^{\text{I\!I}}_{k}(t)=i\psi(t)\zeta_{k}(t).

(113)

This reduces Eq. 112 to,

i\partial_{t}\psi=\hat{H}_{S}\psi-i\lambda^{2}\sum_{j=1}^{K}T_{j}^{\dagger}\chi^{\text{I}}_{j}-i\lambda\sum_{j=1}^{K}T_{j}\chi_{j}^{\text{I\!I}}(t).

(114)

It remains to derive a closed-from equation for Eq. 113. Using the Itô’s formula, we obtain,

i\partial_{t}\chi^{\text{I\!I}}_{k}=(\hat{H}_{S}-d_{k})\chi^{\text{I\!I}}_{k}+i\gamma_{k}\psi(t)\dot{w}_{k}+\lambda\zeta_{k}(t)T_{k}^{\dagger}\chi^{\text{I}}_{k}+\lambda\zeta_{k}(t)T_{k}\chi^{\text{I\!I}}_{k}.

(115)

When the coupling parameter $\lambda$ is sufficiently small, one can add or drop $\mathcal{O}(\lambda)$ terms, which will contribute to an $\mathcal{O}(\lambda^{2})$ error when substituted into Eq. 114. If such an error is acceptable, by collecting equations, we have extended Schrödinger equations,

\left\{\begin{array}[]{l}i\partial_{t}\psi=\displaystyle\hat{H}_{S}\psi-i\lambda\sum_{k=1}^{K}T_{k}^{\dagger}\chi^{\text{I}}_{k}-i\lambda\sum_{k}T_{k}\chi^{\text{I\!I}}_{k},\\ \begin{array}[]{ll}i\partial_{t}\chi^{\text{I}}_{k}=&(\hat{H}_{S}+d^{*}_{k})\chi^{\text{I}}_{k}+i\lambda T_{k}\psi(t),\\ i\partial_{t}\chi^{\text{I\!I}}_{k}=&(\hat{H}_{S}-d_{k})\chi^{\text{I\!I}}_{k}+i\lambda T_{k}^{\dagger}\psi(t)+i\gamma_{k}\psi(t)\dot{w}_{k}.\end{array}\qquad k=1,2,\ldots,K.\end{array}\right.

(116)

Rather than dropping $\mathcal{O}(\lambda)$ terms in Eq. 115, one can continue such a procedure and incorporate the high-order terms. Specifically, noticing the similarity between the $\mathcal{O}(\lambda)$ terms in Eq. 115 with Eq. 113, we define,

\chi_{k}^{\text{I\!I\!I}}=i\zeta_{k}(t)\chi_{k}^{\text{I}},\quad\chi_{k}^{\text{I\!V}}=i\zeta_{k}(t)\chi_{k}^{\text{I\!I}}.

(117)

By repeating the above procedure, one can derive similar equations for these auxiliary wave functions. These embedding steps yield the extended Schrödinger equations (ESE) Eq. 35.

Appendix B The proof of Lemma 9

Proof.

Our proof will mainly target statement (4). The rest of the lemma will become self-evident throughout the proof. In light of the structure of the GQME in Eq. 44, it is enough to consider the case $K=1.$ In this case, $\Gamma$ can be viewed as a $5\times 5$ block matrix. We first write $V_{1}$ in a block matrix form,

V_{1}=\left(\begin{array}[]{cc}0&0\\ R_{1}&0\end{array}\right),

where $R_{1}$ is a $3\times 3$ block matrix with diagonals $\sqrt{2\nu_{1}}I_{S},\sqrt{2\nu_{1}}I_{S}$ and $2\sqrt{2\nu_{1}}I_{S}$ . With direct calculations, we can show that the last term in Eq. 44 can be written as,

V_{1}\Gamma V_{1}^{\dagger}=\left(\begin{array}[]{cc}0&0\\ 0&R_{1}\Gamma_{0:2,0:2}R_{1}^{\dagger}\end{array}\right).

Here $\Gamma_{0:2,0:2}$ refers to the first $3\times 3$ sub-matrix of $\Gamma$ .

We first look at the scenario when $\Gamma(0)$ is block diagonal. The zero blocks in $V_{1}\Gamma V_{1}^{\dagger}$ , along with the observation that $H_{0}$ is block diagonal, imply that the off-diagonals do not change. For the diagonal blocks of $\Gamma(t)$ , we first have,

	$\displaystyle\Gamma_{0,0}(t)=$	$\displaystyle U_{S}(t)\rho_{S}(0)U_{S}(t)^{\dagger},$		(118)
	$\displaystyle\Gamma_{1,1}(t)=$	$\displaystyle\exp(-it(-d_{1}+d_{1}^{*}))U_{S}(t)\Gamma_{1,1}(0)U_{S}(t)^{\dagger}.$		(118)

As a result, these two blocks have norms given respectively by,

\|\Gamma_{0,0}(t)\|=\|\Gamma_{0,0}(0)\|,\quad\|\Gamma_{1,1}(t)\|=\|\Gamma_{1,1}(0)\|e^{-2\nu_{1}t}.

The next three block diagonals will pick up non-homogeneous terms. For instance, we have,

\partial_{t}\Gamma_{2,2}=-i[H_{s}-d_{1},\Gamma_{2,2}]+2\nu_{1}\Gamma_{0,0}.

Using the variation-of-constant formula, we have,

\Gamma_{2,2}(t)=e^{-2\nu_{1}t}U_{S}(t)\Gamma_{2,2}(0)U_{S}(t)^{\dagger}+2\nu_{1}\int_{0}^{t}e^{-2\nu_{1}\tau}U_{S}(\tau)\Gamma_{0,0}(t-\tau)U_{S}(\tau)^{\dagger}d\tau.

(119)

Also by noticing that $\|\Gamma_{1,1}(t)\|$ is constant in time, the diagonal block $\Gamma_{3,3}$ can be bounded directly as,

\|\Gamma_{2,2}(t)\|\leq\|\Gamma_{2,2}(0)\|\exp(-2\nu_{1}t)+\|\Gamma_{0,0}(0)\|(1-e^{-2\nu_{1}t}).

(120)

The right hand side remains bounded for all time. Similarly, the next diagonal block can be expressed as,

	$\displaystyle\Gamma_{3,3}(t)$	$\displaystyle=\exp(-4\nu_{1}t)U_{S}(t)\Gamma_{3,3}(0)U_{S}(t)^{\dagger}$		(121)
		$\displaystyle+2\nu_{1}\int_{0}^{t}\exp(-4\nu_{1}t)U_{S}(\tau)\Gamma_{1,1}(t-\tau)U_{S}(\tau)^{\dagger}d\tau.$		(121)

Notice that $\|\Gamma_{1,1}(t)\|$ is proportional to $\exp(-2\nu_{1}t)$ . Essentially, what leads to the boundedness of the solution is the fact that there is no secular term, implying that $\|\Gamma_{3,3}(t)\|$ follows a similar bound as $\|\Gamma_{2,2}(t)\|_{2}$ . The estimate of $\|\Gamma_{4,4}(t)\|$ follows the same steps.

We now turn to the off-diagonal blocks. By direct calculations, we have, for $j=0,1$ , $k=1,2,3,4$ , and $k>j$

\partial_{t}\Gamma_{j,k}=-i\big{(}H_{S}\Gamma_{j,k}-\Gamma_{j,k}(H_{s}+d_{1})\big{)},

which yields,

\Gamma_{j,k}(t)=\exp(-id_{1}^{*}t)U_{S}(t)\Gamma_{j,k}(0)U_{S}^{\dagger}(t)\Longrightarrow\|\Gamma_{j,k}(t)\|=\|\Gamma_{j,k}(0)\|\exp(-\nu_{1}t).

(122)

For the remaining off-diagonal entries, we will check $\Gamma_{2,3}$ as an example. It follows the equation,

\partial_{t}\Gamma_{2,3}=-i\big{(}(H_{s}-d_{1})\Gamma_{2,3}-\Gamma_{2,3}(H_{s}+d_{1}+d_{1}^{*})\big{)}+2\nu_{1}\Gamma_{0,1}.

This implies that,

\|\Gamma_{3,4}(t)\|\leq\|\Gamma_{3,4}(0)\|e^{-3\nu_{1}t}+\|\Gamma_{1,2}(0)\|2e^{-2\nu_{1}t}(1-e^{-\nu_{1}t}).

We also see from these calculations that these off-diagonal blocks will become zero if the initial matrix $\Gamma(0)$ is block diagonal.

By examining the block entries of $\Gamma(t)$ , we have shown the boundedness of the solution stated in Lemma 9 for all time. ∎

Appendix C The Proof of Theorem 10

Proof.

We will prove the asymptotic bound using an expansion of the Eq. 27. More specifically, we write the total density matrix in terms of powers of $\lambda,$

\rho(t)=\rho^{(0)}(t)+\lambda\rho^{(1)}(t)+\lambda^{2}\rho^{(2)}(t)+\mathcal{O}(\lambda^{3}).

(123)

By taking a partial trace over the bath space, we obtain a similar expansion for $\rho_{S}:$

\rho_{S}(t)=\rho_{S}^{(0)}(t)+\lambda\rho_{S}^{(1)}(t)+\lambda^{2}\rho_{S}^{(2)}(t)+\mathcal{O}(\lambda^{3}).

(124)

By inserting Eq. 123 into Eq. 27 and separate terms of different order, we arrive at

$\displaystyle i\partial_{t}\rho^{(0)}$	$\displaystyle=[H_{S}\otimes I_{B}+I_{S}\otimes H_{B},\rho^{(0)}],\;\;\rho^{(0)}(0)=\rho(0),$	(125)
$\displaystyle i\partial_{t}\rho^{(1)}$	$\displaystyle=[H_{S}\otimes I_{B}+I_{S}\otimes H_{B},\rho^{(1)}]+\sum_{\alpha=1}^{M}[S_{\alpha}\otimes B_{\alpha},\rho^{(0)}],\;\;\rho^{(1)}(0)=0,$
$\displaystyle i\partial_{t}\rho^{(2)}$	$\displaystyle=[H_{S}\otimes I_{B}+I_{S}\otimes H_{B},\rho^{(2)}]+\sum_{\alpha=1}^{M}[S_{\alpha}\otimes B_{\alpha},\rho^{(1)}],\;\;\rho^{(2)}(0)=0.$

Within this expansion, the dynamics of $\rho^{(0)}$ contains no coupling. Let

U(t)=U_{S}(t)\otimes U_{B}(t),\quad U_{S}=\exp\left(-itH_{S}\right),\;U_{B}(t)=\exp\left(-itH_{B}\right),

be the unitary operators. Then we have,

\rho^{(0)}(t)=U(t)\rho^{(0)}(0)U(t)^{\dagger}.

(126)

Since all the operators on the right hand side are in tensor product forms, $\rho^{(0)}(t)$ remains a tensor product:

\rho^{(0)}(t)=\rho_{S}^{(0)}(t)\otimes\rho_{B}.

(127)

Here $\rho_{S}^{(0)}(t)=U_{S}(t)\rho_{S}(0)U_{S}(t)^{\dagger}.$ Meanwhile, the matrix $\rho_{B}$ stays because it commutes with $H_{B}.$

The term $\rho^{(1)}$ can be expressed using the variation-of-constant formula,

	$\displaystyle\rho^{(1)}(t)=$	$\displaystyle-i\sum_{\alpha}\int_{0}^{t}U(t-t^{\prime})S_{\alpha}\otimes B_{\alpha}\rho^{(0)}(t^{\prime})U(t-t^{\prime})^{\dagger}dt^{\prime}$
		$\displaystyle+i\sum_{\alpha}\int_{0}^{t}U(t-t^{\prime})\rho^{(0)}(t^{\prime})S_{\alpha}\otimes B_{\alpha}U(t-t^{\prime})^{\dagger}dt^{\prime}.$

In light of Eq. 127, we can make the same observation that $\rho^{(1)}(t)$ consists of terms that are tensor products. By following standard notations [BP02], i.e.,

S_{\alpha}(t)=U_{S}(t)^{\dagger}S_{\alpha}U_{S}(t),\quad B_{\alpha}(t)=U_{B}(t)^{\dagger}B_{\alpha}U_{B}(t),

(128)

we can simplify $\rho^{(1)}(t)$ as follows,

	$\displaystyle\rho^{(1)}(t)=$	$\displaystyle-i\sum_{\alpha}\int_{0}^{t}S_{\alpha}(t^{\prime}-t)\rho_{S}^{(0)}(t)\otimes B_{\alpha}(t^{\prime}-t)\rho_{B}dt^{\prime}$
		$\displaystyle+i\sum_{\alpha}\int_{0}^{t}\rho_{S}^{(0)}(t)S_{\alpha}(t^{\prime}-t)\otimes B_{\alpha}(t^{\prime}-t)\rho_{B}dt^{\prime}.$

Since $\rho_{B}$ commutes with $H_{B}$ , it commutes with $U_{B}.$ Therefore,

\tr(B_{\alpha}(t^{\prime}-t)\rho_{B})=\tr(B_{\alpha}\rho_{B})=0,

which shows that

\tr_{B}\left(\rho^{(1)}(t)\right)=0.

Therefore, $\rho_{S}^{(1)}(t)=0$ . The correction to $\rho_{S}$ comes from $\rho_{S}^{(2)}(t),$ which is similarly expressed as,

	$\displaystyle\rho^{(2)}(t)\!=$	$\displaystyle-i\sum_{\alpha}\int_{0}^{t}U(t-t^{\prime})S_{\alpha}\otimes B_{\alpha}\rho^{(1)}(t^{\prime})U(t-t^{\prime})^{\dagger}dt^{\prime}$
		$\displaystyle+i\sum_{\alpha}\int_{0}^{t}U(t-t^{\prime})\rho^{(1)}(t^{\prime})S_{\alpha}\otimes B_{\alpha}U(t-t^{\prime})^{\dagger}dt^{\prime},$
	$\displaystyle=$	$\displaystyle\!-\!\sum_{\alpha}\!\sum_{\beta}\!\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\alpha}(t^{\prime}-t)S_{\beta}(\tau-t)\rho_{S}^{(0)}(t)\!\otimes\!B_{\alpha}(t^{\prime}-t)B_{\beta}(\tau-t)\rho_{B}d\tau dt^{\prime},$
		$\displaystyle\!+\!\sum_{\alpha}\!\sum_{\beta}\!\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\beta}(\tau-t)\rho_{S}^{(0)}(t)S_{\alpha}(t^{\prime}-t)\!\otimes\!B_{\beta}(\tau-t)\rho_{B}B_{\alpha}(t^{\prime}-t)d\tau dt^{\prime},$
		$\displaystyle\!+\!\sum_{\alpha}\!\sum_{\beta}\!\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\alpha}(t^{\prime}-t)\rho_{S}^{(0)}(t)S_{\beta}(\tau-t)\!\otimes\!B_{\alpha}(t^{\prime}-t)\rho_{B}B_{\beta}(\tau-t)d\tau dt^{\prime},$
		$\displaystyle\!-\!\sum_{\alpha}\!\sum_{\beta}\!\int_{0}^{t}\int_{0}^{t^{\prime}}\rho_{S}^{(0)}(t)S_{\beta}(\tau-t)S_{\alpha}(t^{\prime}-t)\!\otimes\!\rho_{B}B_{\beta}(\tau-t)B_{\alpha}(t^{\prime}-t)d\tau dt^{\prime}.$

Invoking the bath correlation function,

C_{\alpha,\beta}(t)=\tr(B_{\alpha}(t)B_{\beta}\rho_{B}),

(129)

we arrive at an expansion of $\rho_{S}(t)$ up to ${\mathcal{O}\left(\lambda^{2}\right)}$ ,

$\displaystyle\rho_{S}(t)=\rho_{S}^{(0)}(t)$	$\displaystyle-\lambda^{2}\sum_{\alpha}\sum_{\beta}\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\alpha}(t^{\prime}-t)S_{\beta}(\tau-t)\rho_{S}^{(0)}(t)C_{\alpha,\beta}(t^{\prime}-\tau)d\tau dt^{\prime}$	(130)
	$\displaystyle+\lambda^{2}\sum_{\alpha}\sum_{\beta}\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\beta}(\tau-t)\rho_{S}^{(0)}(t)S_{\alpha}(t^{\prime}-t)C_{\alpha,\beta}(t^{\prime}-\tau)d\tau dt^{\prime}$
	$\displaystyle+\lambda^{2}\sum_{\alpha}\sum_{\beta}\int_{0}^{t}\int_{0}^{t^{\prime}}S_{\alpha}(t^{\prime}-t)\rho_{S}^{(0)}(t)S_{\beta}(\tau-t)C_{\beta,\alpha}(t^{\prime}-\tau)^{*}d\tau dt^{\prime}$
	$\displaystyle-\lambda^{2}\sum_{\alpha}\sum_{\beta}\int_{0}^{t}\int_{0}^{t^{\prime}}\rho_{S}^{(0)}(t)S_{\beta}(\tau-t)S_{\alpha}(t^{\prime}-t)C_{\alpha,\beta}(t^{\prime}-\tau)^{*}d\tau dt^{\prime}$
	$\displaystyle+\mathcal{O}(\lambda^{3}).$

Here we have used the property of the bath correlation function: $C_{\alpha,\beta}(t)=C_{\beta,\alpha}(-t)^{*}.$ Now we incorporate the function form of the bath correlation function in Eq. 5. We find that,

$\displaystyle\rho_{S}(t)=\rho_{S}^{(0)}(t)$	$\displaystyle-\lambda^{2}\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}T_{k}(t^{\prime}-t)T_{k}(\tau-t)\rho_{S}^{(0)}(t)e^{-i(t^{\prime}-\tau)d_{k}t}d\tau dt^{\prime}$	(131)
	$\displaystyle+\lambda^{2}\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}T_{k}(\tau-t)\rho_{S}^{(0)}(t)T_{k}(t^{\prime}-t)e^{-i(t^{\prime}-\tau)d_{k}t}d\tau dt^{\prime}$
	$\displaystyle+\lambda^{2}\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}T_{k}(t^{\prime}-t)\rho_{S}^{(0)}(t)T_{k}(\tau-t)e^{i(t^{\prime}-\tau)d_{k}t}d\tau dt^{\prime}$
	$\displaystyle-\lambda^{2}\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}\rho_{S}^{(0)}(t)T_{k}(\tau-t)T_{k}(t^{\prime}-t)e^{i(t^{\prime}-\tau)d_{k}t}d\tau dt^{\prime}$
	$\displaystyle+\mathcal{O}(\lambda^{3}).$

Now we show that the GQME in Eq. 44 has an asymptotic expansion that is consistent with Eq. 130. Expanding $\Gamma$ as,

\Gamma(t)=\Gamma^{(0)}(t)+\lambda\Gamma^{(1)}(t)+\lambda^{2}\Gamma^{(2)}(t)+\mathcal{O}(\lambda^{3}),

and substituting it into Eq. 44, one gets,

$\displaystyle\Gamma^{(0)}(t)$	$\displaystyle=\exp\left(t\mathcal{L}_{0}\right)\Gamma(0),$	(132)
$\displaystyle\Gamma^{(1)}(t)$	$\displaystyle=-i\int_{0}^{t}\exp\left((t-t^{\prime})\mathcal{L}_{0}\right)[H_{1},\Gamma^{(0)}(t^{\prime})]dt^{\prime},$
$\displaystyle\Gamma^{(2)}(t)$	$\displaystyle=-i\int_{0}^{t}\exp\left((t-t^{\prime})\mathcal{L}_{0}\right)[H_{1},\Gamma^{(1)}(t^{\prime})]dt^{\prime}.$

The leading term $\Gamma^{(0)}(t)$ has been shown to be a block diagonal matrix in the previous section. The first diagonal block is precisely $\rho_{S}^{(0)}(t)$ , which is consistent with the ${\mathcal{O}\left(1\right)}$ term in Eq. 130. To examine $\Gamma^{(1)}(t)$ , we first notice that the commutator in the integral has the following structure,

\Xi^{(0)}(t^{\prime})=[H_{1},\Gamma^{(0)}(t^{\prime})]=i\left[\begin{array}[]{cccccc}0&\Xi_{0,1}^{(0)}(t^{\prime})&\Xi_{0,2}^{(0)}(t^{\prime})&0&0&\cdots\\ \Xi_{1,0}^{(0)}(t^{\prime})&0&0&0&0&\cdots\\ \Xi_{2,1}^{(0)}(t^{\prime})&0&0&\Xi_{3,4}^{(0)}(t^{\prime})&\Xi_{3,5}^{(0)}(t^{\prime})&\cdots\\ 0&0&\Xi_{4,3}^{(0)}(t^{\prime})&0&0&\cdots\\ 0&0&\Xi_{5,3}^{(0)}(t^{\prime})&0&0&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right].

(133)

Here we highlighted the leading $5\times 5$ submatrix and the zero blocks within it. This is enough for the purpose of the proof.

By inspecting the solutions that correspond to $\exp\left(t\mathcal{L}_{0}\right)$ , we find that the first diagonal block of $\int_{0}^{t}\exp\left((t-t^{\prime})\mathcal{L}_{0}\right)\Xi^{(0)}(t^{\prime})$ is zero. Therefore, $\Gamma^{(1)}$ has no contribution to the density matrix $\rho_{S}$ , i.e., there is no ${\mathcal{O}\left(\lambda\right)}$ term. This is consistent with Eq. 130.

To proceed further, we have to identify the nonzero blocks in Eq. 133. With direct calculations, we have,

	$\displaystyle\Xi_{0,4k-3}^{(0)}(t)=$	$\displaystyle-\Gamma_{1,1}^{(0)}(t)T_{j}=-\rho_{S}^{(0)}(t)T_{j},$
	$\displaystyle\Xi_{0,4k-2}^{(0)}(t)=$	$\displaystyle T_{k}\Gamma_{4k-2,4k-2}^{(0)}(t).$

From Eq. 132, we can extract the equation,

i\partial_{t}\Gamma_{0,4k-3}^{(1)}=H_{S}\Gamma_{0,4k-3}^{(1)}-\Gamma_{0,4k-3}^{(1)}(H_{S}+d_{k}).

Combining the two equations above, we obtain,

	$\displaystyle\Gamma_{0,4k-3}^{(1)}(t^{\prime})$	$\displaystyle=\int_{0}^{t^{\prime}}U_{S}(t^{\prime}-t^{\prime\prime})\rho_{S}^{(0)}(t^{\prime\prime})T_{k}U_{S}(t^{\prime}-t^{\prime\prime})^{\dagger}e^{i(t^{\prime}-t^{\prime\prime})d_{k}}dt^{\prime\prime}$
		$\displaystyle=\int_{0}^{t^{\prime}}\rho_{S}^{(0)}(t^{\prime})T_{k}(t^{\prime\prime}-t^{\prime})e^{i(t^{\prime}-t^{\prime\prime})d_{k}}dt^{\prime\prime}.$

Again using the solution properties associated with the superoperator $\exp\left(t\mathcal{L}_{0}\right)$ , we have that the first block of $\Gamma^{(1)}$ is given by,

		$\displaystyle-i\int_{0}^{t}U_{S}(t-t^{\prime})\Xi^{(1)}(t^{\prime})U_{S}(t-t^{\prime})^{\dagger}dt^{\prime}$
	$\displaystyle=$	$\displaystyle\sum_{k}\int_{0}^{t}U_{S}(t-t^{\prime})\big{(}[\Gamma_{0,4k-3}(t^{\prime}),T_{k}]+[\Gamma_{0,4k-2}(t^{\prime}),T_{k}^{\dagger}]\big{)}U_{S}(t-t^{\prime})^{\dagger}dt^{\prime}.$

Similar to Eq. 131, we have also obtained four terms after expanding the commutators. Let us examine the first integral term,

		$\displaystyle\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}U_{S}(t-t^{\prime})\rho_{S}^{(0)}(t^{\prime})T_{k}(t^{\prime\prime}-t^{\prime})T_{k}e^{i(t^{\prime}-t^{\prime\prime})d_{k}}U_{S}(t-t^{\prime})^{\dagger})dt^{\prime\prime}$
		$\displaystyle=\sum_{k}\int_{0}^{t}\int_{0}^{t^{\prime}}\rho_{S}^{(0)}(t^{\prime})T_{k}(t^{\prime\prime}-t)T_{k}(t^{\prime}-t)e^{i(t^{\prime}-t^{\prime\prime})d_{k}}dt^{\prime\prime}dt^{\prime}.$

This is the same as the first last integral in Eq. 131. The rest of the integrals can be similarly verified. ∎

Appendix D The proof of Theorem 11

Proof.

In Lemma 9, we have proved a bound for the case when $\lambda=0.$ Denote the solution by $\Gamma_{0}(t)$ , and the solution operator by $\exp(t\mathcal{L}_{0})$ . Therefore, the GQME in Eq. 44 can be written in a perturbation form,

\partial_{t}\Gamma=\mathcal{L}_{0}\Gamma-i\lambda(H_{1}\Gamma-\Gamma H_{1}^{\dagger}).

(134)

The solution can be recast in an integral form,

\Gamma(t)=\Gamma_{0}(t)-i\lambda\int_{0}^{t}\exp\left((t-\tau)\mathcal{L}_{0}\right)(H_{1}\Gamma-\Gamma H_{1}^{\dagger})d\tau.

(135)

From Lemma 8, the super-operator $\exp(t\mathcal{L}_{0})$ is bounded. Therefore, we have,

\|\Gamma(t)-\Gamma_{0}(t)\|\leq 2\lambda C\|H_{1}\|\int_{0}^{t}\|\Gamma(t)\|dt.

As a result, the bound can be obtained by directly using Gronwall’s inequality.

Since $\Gamma_{0}(t)$ is block diagonal, the above inequality also shows that the off-diagonal blocks of $\Gamma(t)$ are of order $\lambda.$ Finally, the bounds for the trace can be verified from Eq. 118, Eq. 119 and Eq. 121 in the proof of Lemma 9. ∎

Appendix E The proof of Corollary 12

Proof.

From Eq. 134, we may take the trace.

\tr(\Gamma(t))-\tr(\Gamma^{(0)}(t))=-i\lambda\int_{0}^{t}\tr(\Sigma(t,\tau))d\tau,

(136)

where,

\Sigma(t,\tau)=\exp(t\mathcal{L}_{0})\Xi,

with $\Xi$ from Eq. 52. By the definition of $\exp(t\mathcal{L}_{0})$ , $\Sigma(t,\tau)$ is the solution of the equation,

\partial_{t}\Sigma=\mathcal{L}_{0}\Sigma,\quad\Sigma(\tau,\tau)=(H_{1}\Gamma-\Gamma H_{1}^{\dagger}).

Using the property in Eq. 59 of the super-operator $\exp(t\mathcal{L}_{0})$ in Theorem 11, we obtain the bound,

\mid\tr\big(\Sigma(t,\tau)\big{missing})\mid\leq 3K\mid\tr\big(\Sigma_{0,0}(\tau,\tau)\big{missing})\mid+\sum_{k>0}\mid\tr\big(\Sigma_{k,k}(\tau,\tau)\big{missing})\mid.

We now invoke the estimate in Lemma 8. The trace of each diagonal block of $\Sigma(\tau,\tau)$ is bounded by the off-diagonal blocks of $\Gamma(t)$ , which is of order $\lambda$ . Namely, there exists a constant, such that,

\mid\tr\big(\Sigma(\tau,\tau)\big{missing})\mid\leq C\lambda.

Collecting terms, we have,

\mid\tr(\Gamma(t))-\tr(\Gamma^{(0)}(t))\mid\leq Ct\lambda^{2}.

Alternatively, one can start with Eq. 134, and apply the formula again to replace the $\Gamma(t)$ in the integral:

	$\displaystyle\Gamma(t)-\Gamma^{(0)}(t)=$	$\displaystyle-i\lambda\int_{0}^{t}\exp\left((t-\tau)\mathcal{L}_{0}\right)(H_{1}\Gamma^{(0)}(\tau)-\Gamma^{(0)}(\tau)H_{1}^{\dagger})d\tau$
		$\displaystyle-\lambda^{2}\int_{0}^{t}\int_{0}^{\tau}\exp\left((t-\tau)\mathcal{L}_{0}\right)(H_{1}\Sigma(\tau,s)-\Sigma(\tau,s)H_{1}^{\dagger})dsd\tau+\mathcal{O}(\lambda^{3}).$

Here $\Sigma(\tau,s)=\exp\left((\tau-s)\mathcal{L}_{0}\right)\Xi_{0}$ , with

\Xi^{(0)}=H_{1}\Gamma_{0}(s)-\Gamma_{0}(s)H_{1}^{\dagger}.

(137)

For the $\mathcal{O}(\lambda)$ term, we notice that $\Gamma{(0)}$ is block diagonal, and so in light of Lemma 8, the matrix $\Xi^{(0)}$ has zero trace. This implies that,

\tr\big(\Gamma(t))-\Gamma^{(0)}(t)\big{)}=\mathcal{O}(\lambda^{2}).

For the $\mathcal{O}(\lambda^{2})$ term, from Eq. 55, we have, the trace of the first block of $H_{1}\Sigma(\tau,s)-\Sigma(\tau,s)H_{1}^{\dagger}$ is given by,

\sum_{k=1}^{K}\tr\Big((\Sigma_{0,4k-3}-\Sigma_{4k-2,0})T_{k}\Big{missing})+\sum_{k=1}^{K}\tr\Big(T_{k}^{\dagger}(\Sigma_{0,4k-2}-\Sigma_{4k-3,0})\Big{missing}).

(138)

Meanwhile, from the proof of Lemma 9, we see that the superoperator does not change the off-diagonal blocks in the first row and column. Thus, $\Sigma_{0,j}(\tau,s)=\Sigma_{0,j}(s,s)=\Xi_{0,j}^{(0)}(s)$ .

With direct matrix multiplications, we find that,

\Xi_{4k-3,0}=-T_{k}\Gamma_{4k-3,4k-3}-\Gamma_{0,0}T_{k}^{\dagger},\quad\Xi_{4k-2,0}=T_{k}\Gamma_{4k-2,4k-2}.

From the proof of Lemma 9, we also have $\Gamma_{4k-3,4k-3}(t)=0.$ Combining these steps, we find that the trace in Eq. 138 is zero. Therefore, we have

\tr(\rho_{S}(t))=1+\mathcal{O}(\lambda^{3}).

∎

Appendix F The proof of Lemma 16

Proof.

We use an intermediate superoperator $\mathcal{I}+\delta\mathcal{K}$ . Assume the Hilbert space $\mathcal{K}$ is acting on has dimension $N$ . Consider a Hilbert space of arbitrary dimension $N^{\prime}$ . For any operator $Q$ acting on $\mathbb{C}^{N}\otimes\mathbb{C}^{N^{\prime}}$ with $\norm{Q}_{1}=1$ , we have

		$\displaystyle\norm{(\mathcal{M}_{\delta}\otimes\mathcal{I}_{N^{\prime}}-(\mathcal{I}_{N}+\delta\mathcal{K})\otimes\mathcal{I}_{N^{\prime}})(Q)}_{1}$		(139)
		$\displaystyle=\norm{\sum_{j=0}^{m}(A_{j}\otimes I)Q(A_{j}\otimes I)^{{\dagger}}-(Q+\delta(\mathcal{K}\otimes\mathcal{I}_{N^{\prime}})(Q)}_{1}$
		$\displaystyle=\norm{\delta^{2}(H\otimes I)Q(H^{{\dagger}}\otimes I)}_{1}$
		$\displaystyle\leq\norm{H\otimes I}^{2}$
		$\displaystyle\leq(\delta\Lambda)^{2}.$

Now, we have

\displaystyle\norm{(\mathcal{M}_{\delta}-(\mathcal{I}_{N}+\delta\mathcal{K})}_{\diamond}\leq(\delta\Lambda)^{2}.

(140)

To bound the distance between $\mathcal{I}+\delta\mathcal{K}$ and $e^{\delta\mathcal{K}}$ , we assume $0\leq\delta\norm{\mathcal{K}}_{\diamond}\leq 1$ . Consider any $X$ such that $\norm{X}_{1}\leq 1$ , we have

	$\displaystyle\norm{(e^{\delta\mathcal{K}}-(\mathcal{I}+\delta\mathcal{K}))(X)}_{1}$	$\displaystyle=\norm{\sum_{s=2}^{\infty}\frac{\delta^{s}}{s!}\mathcal{K}^{s}(X)}_{1}\leq\sum_{s=2}^{\infty}\frac{\delta^{s}}{s!}\norm{\mathcal{K}^{s}(X)}_{1}$		(141)
		$\displaystyle\leq\sum_{s=2}^{\infty}\frac{\delta^{s}}{s!}\norm{\mathcal{K}(X)}_{1}^{s}\leq(\delta\norm{\mathcal{K}(X)}_{1})^{2}\leq(\delta\norm{\mathcal{K}}_{1})^{2},$		(142)

where the penultimate inequality follows from the fact that $e^{z}-(1+z)\leq z^{2}$ when $0\leq z\leq 1$ .

Now, we extend this bound to the diamond norm. Note that, for two Hilbert spaces $\mathbb{C}^{N}$ and $\mathbb{C}^{N^{\prime}}$ ,

\displaystyle(e^{\delta\mathcal{K}}-(\mathcal{I}_{N}+\delta\mathcal{K}))\otimes\mathcal{I}_{N^{\prime}}=e^{\delta(\mathcal{K}\otimes\mathcal{I}_{N^{\prime}})}-(\mathcal{I}_{N\times N^{\prime}}+\delta(\mathcal{K}\otimes\mathcal{I}_{N^{\prime}})).

(143)

When $N=N^{\prime}$ , we have that $\norm{\mathcal{K}\otimes\mathcal{I}_{N^{\prime}}}_{1}=\norm{\mathcal{K}}_{\diamond}$ . This implies that

$\displaystyle\norm{(e^{\delta\mathcal{K}}-(\mathcal{I}_{N^{\prime}}+\delta\mathcal{K})}_{\diamond}$	$\displaystyle=\norm{\left(e^{\delta\mathcal{K}}-(\mathcal{I}_{N}+\delta\mathcal{K})\right)\otimes\mathcal{I}_{N^{\prime}}}_{1}$	(144)
	$\displaystyle=\norm{e^{\delta(\mathcal{K}\otimes\mathcal{I}_{N^{\prime}})}-(\mathcal{I}_{N\times N^{\prime}}+\delta(\mathcal{K}\otimes\mathcal{I}_{N^{\prime}}))}_{1}$	(145)
	$\displaystyle\leq(\delta\norm{\mathcal{K}\otimes\mathcal{I}_{N^{\prime}}}_{1})^{2}$	(146)
	$\displaystyle\leq(\delta\norm{\mathcal{K}}_{\diamond})^{2}.$	(147)

To see the relationship between $\norm{\mathcal{L}}_{\diamond}$ and $\norm{\mathcal{K}}_{1}$ , first observe that $\norm{\mathcal{K}}_{1}\leq 2\Lambda$ . Then using the fact the $\norm{M\otimes I}=\norm{M}$ for all $M$ , it follows that $\norm{\mathcal{K}}_{\diamond}\leq 2\Lambda$ . Together with Eq. 147, we have

\displaystyle\norm{(e^{\delta\mathcal{K}}-(\mathcal{I}_{N^{\prime}}+\delta\mathcal{K})}_{\diamond}\leq(2\delta\Lambda)^{2}.

(148)

This result in the lemma now directly follows from Eqs. 140 and 148. ∎

$\displaystyle Q\ket{\psi}$	$\displaystyle=(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)W^{{\dagger}}P_{0}W(\ket{\widehat{0}}\ket{\widehat{\mu}}\otimes\ket{\psi})=\sin\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)W^{{\dagger}}\ket*{\widehat{\phi}}$	(82)
	$\displaystyle=\sin\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)\left(\sin\theta\ket{\widehat{\psi}}+\cos\theta\ket{\widehat{\psi}^{\bot}}\right)$
	$\displaystyle=\sin^{2}\theta\ket{\psi}+\sin\theta\cos\theta(\bra{\widehat{0}}\bra{\widehat{\mu}}\otimes I)\ket*{\widehat{\psi}^{\bot}},$