Hardware-efficient quantum algorithm for the simulation of open-system dynamics and thermalisation

TL;DR

This paper introduces a quantum algorithm that efficiently simulates open-system dynamics and thermalisation using a minimal number of qubits by reproducing reservoir correlation functions, enabling practical quantum simulations of complex environments.

Contribution

The authors propose a novel method to simulate open-system dynamics with fewer qubits by reproducing reservoir correlation functions, reducing resource requirements for quantum simulations.

Findings

Efficient simulation of open-system dynamics with fewer qubits.

Reproduction of up to n-time correlation functions for master equation simulation.

Environment size can be smaller than the system in the algorithm.

Abstract

The quantum open-system simulation is an important category of quantum simulation. By simulating the thermalisation process at the zero temperature, we can solve the ground-state problem of quantum systems. To realise the open-system evolution on the quantum computer, we need to encode the environment using qubits. However, usually the environment is much larger than the system, i.e. numerous qubits are required if the environment is directly encoded. In this paper, we propose a way to simulate open-system dynamics by reproducing reservoir correlation functions using a minimised Hilbert space. In this way, we only need a small number of qubits to represent the environment. To simulate the $n$-th-order expansion of the time-convolutionless master equation by reproducing up to $n$-time correlation functions, the number of qubits representing the environment is $\sim \lfloor \frac{n}{2} \rfloor \log_2(N_\omega N_\beta)$. Here, $N_\omega$ is the number of frequencies in the discretised environment spectrum, and $N_\beta$ is the number of terms in the system-environment interaction. By reproducing two-time correlation functions, i.e. taking $n = 2$, we can simulate the Markovian quantum master equation. In our algorithm, the environment on the quantum computer could be even smaller than the system.

Equations

$$\frac{\partial }{\partial t}\rho (t)=-i\alpha [H_{\mathrm{I}}(t),\rho (t)]\equiv \alpha \mathcal{L}(t)\rho (t).$$

$$\frac{\partial }{\partial t}\mathcal{P}\rho (t)=\mathcal{K}(t)\mathcal{P}\rho (t)$$

$${\mathcal{K}}_4(t)$$

$${\mathbf{\kappa }}_n(t,t_1,\dots ,t_{n-1})=\mathcal{P}\mathcal{L}(t)\mathcal{L}(t_1)\cdots \mathcal{L}(t_{n-1})\mathcal{P}.$$

$${\mathbf{\kappa }}_n(t,\dots ,t_{n-1})\rho$$

$$\displaystyle\mathrm{Tr}\big{[}\widetilde{\mathcal{B}}_{\beta}(t,\nu)\cdots\widetilde{\mathcal{B}}_{\beta_{m-1}}(t_{m-1},\nu_{m-1})\widetilde{\rho}_{\rm E}\big{]}$$

$$\frac{d{\rho }_{\mathrm{S}}}{dt}=-{\alpha }^2{\int }_0^{t}ds{\mathrm{Tr}}_{\mathrm{E}}[H_{\mathrm{I}}(t),[H_{\mathrm{I}}(t-s),{\rho }_{\mathrm{S}}\otimes {\rho }_{\mathrm{E}}]].$$

$$\displaystyle\langle B_{\beta}(t)B_{\beta^{\prime}}(t-s)\rangle\equiv\mathrm{Tr}\big{[}B_{\beta}(t)B_{\beta^{\prime}}(t-s)\rho_{\rm E}\big{]}.$$

$$⟨B_{\beta }(t)B_{{\beta }^{\prime }}(t-s)⟩=\sum _{\omega }{e}^{-i\omega s}{\gamma }_{\beta ,{\beta }^{\prime }}(\omega ).$$

$${\tilde{H}}_{\mathrm{E}}$$

$${\tilde{B}}_{\beta }$$

$${\tilde{B}}_{\beta }(\omega )=\sum _l{g}_{\beta ,l}(\omega ){\sigma }_l(\omega )+g_{\beta ,l}^{\ast }(-\omega ){\sigma }_l^{†}(-\omega ).$$

$$\displaystyle\mathrm{Tr}\big{[}\mathcal{B}_{\beta}(t,\nu)\cdots\mathcal{B}_{\beta_{m-1}}(t_{m-1},\nu_{m-1})\rho_{\rm E}\big{]}$$

$$\mathrm{Tr}[\text{(some }B\text{ operators)}{\rho }_{\mathrm{E}}\text{(some }B\text{ operators)}],$$

$$\mathrm{Tr}[\text{(some }b\text{ operators)}|\psi ⟩⟨\psi |\text{(some }b\text{ operators)}].$$

$$\displaystyle\mathrm{Tr}\big{[}\mathcal{B}_{\beta}(\omega,\nu)\cdots\mathcal{B}_{\beta_{m-1}}(\omega_{m-1},\nu_{m-1})\rho_{\rm E}\big{]}$$

$$|{\varphi }_{\Omega ,m}(\beta ,\omega ,{\beta }_1,{\omega }_1,\dots ,{\beta }_{m-1},{\omega }_{m-1})⟩$$

$$\displaystyle\mathrm{Tr}\big{[}\mathcal{B}_{\beta}(\omega,\nu)\cdots\mathcal{B}_{\beta_{m-1}}(\omega_{m-1},\nu_{m-1})\rho_{\rm E}\big{]}$$

$$⟨\tilde{\psi }|{\tilde{b}}_m\cdots {\tilde{b}}_2{\tilde{b}}_1|\tilde{\psi }⟩=⟨\psi |b_m\cdots b_2{b}_1|\psi ⟩$$

$${\tilde{H}}_{\mathrm{E}}$$

$${\tilde{B}}_{\beta }$$

$$\mathrm{Tr}[{\tilde{B}}_{\beta }(t){\tilde{B}}_{{\beta }_3}(t_3){\tilde{\rho }}_{\mathrm{E}}{\tilde{B}}_{{\beta }_2}(t_2){\tilde{B}}_{{\beta }_1}(t_1)]$$

$$\mathrm{Tr}[{\tilde{B}}_{\beta }(t){\tilde{\rho }}_{\mathrm{E}}{\tilde{B}}_{{\beta }_1}(t_1)]\times \mathrm{Tr}[{\tilde{B}}_{{\beta }_3}(t_3){\tilde{\rho }}_{\mathrm{E}}{\tilde{B}}_{{\beta }_2}(t_2)]$$

$$\displaystyle\mathrm{Tr}\big{[}\widetilde{\mathcal{B}}_{\beta}(t,\nu)\cdots\widetilde{\mathcal{B}}_{\beta_{m-1}}(t_{m-1},\nu_{m-1})\widetilde{\rho}_{\rm E}\big{]}$$

$${\tilde{B}}_{\beta }=\sum _{x,l}{g}_{\beta ,l}(x){\sigma }_l(x)+h.c.,$$

$$g_{\beta ,l}(x)=\frac{1}{\sqrt{N_{\omega }}}\sum _{\omega }{e}^{-i\frac{2\pi kx}{N_{\omega }}}{g}_{\beta ,l}(\omega ).$$

$$⟨{\tilde{B}}_{\beta }(t){\tilde{B}}_{{\beta }^{\prime }}(t-s)⟩=⟨\mathrm{v}|{\tilde{B}}_{\beta }{e}^{-i{\tilde{H}}_{\mathrm{E}}s}{\tilde{B}}_{{\beta }^{\prime }}|\mathrm{v}⟩$$

$$\tilde{B}=\sum _{\omega }\left[g(\omega )|\omega ⟩⟨\mathrm{v}|+g^{\ast }(\omega )|\mathrm{v}⟩⟨\omega |\right],$$

$$\Delta$$

Full text

I Introduction

The idea of quantum computation is motivated by quantum simulation. According to R. Feynman, “the physical world is quantum mechanical, and therefore the proper problem is the simulation of quantum physics” Feynman1982 . The physical world is not only quantum but also open. Many vital phenomena are attributed to the open-system dynamics, e.g. thermalisation BreuerPetruccione ; Vega2017 . Systems are influenced by their environments through external interactions. Therefore, by simulating the composite system, including the system and the environment, we can study an open system on a quantum computer Lloyd1996 ; Terhal2000 ; Wang2011 . However, the simulation of the environment is usually inefficient when the environment is big compared to the system. It is also a waste of resources. In many circumstances, we are only interested in the system, not the environment. The simulation of the environment using the most of computational resources may not give us any new knowledge, for instance, when the environment is modeled as exactly solvable boson bath or spin bath. The dynamics of the system is determined by reservoir correlation functions. For example, in the thermalisation, transition rates between eigenstates are determined by two-time reservoir correlation functions BreuerPetruccione . Therefore, reproducing reservoir correlation functions is sufficient, and the full simulation of the environment is unnecessary.

An application of quantum computation is to compute the ground-state energy, which is an important problem in material science and chemistry Abrams1999 ; AspuruGuzik2005 ; Wecker2014 ; Bauer2016 . Given an initial state with a finite probability in the ground state, we can use the quantum phase estimation algorithm to obtain the ground-state energy Abrams1999 ; AspuruGuzik2005 . However, we do not have a universal algorithm that can prepare such an initial state Verstraete2009 ; Farhi2001 . Solving the ground-state problem for a general Hamiltonian is likely to be intractable even in quantum computation Kitaev2002 ; Aharonov2002 . A related problem is preparing or sampling thermal states of a quantum system Poulin2009 ; Bilgin2010 ; Temme2011 ; Riera2012 ; Yung2012 ; Motta2019 , and the ground state is the thermal state at the zero temperature. If we only focus on systems in the real world, most of them reach the thermal state as a result of the open-system dynamics. Therefore, for such real-world systems, simulating the open-system dynamics is an efficient way to prepare thermal states, including the ground state. Although we have quantum algorithms that can implement semi-group dynamics (unitary or non-unitary) Lloyd1996 ; Berry2007 ; Wiebe2010 ; Berry2015 ; Campbell2019 ; Bacon2001 ; Kliesch2011 ; Sweke2015 ; Candia2015 ; Sweke2016 ; Childs2017 ; Chenu2017 , they cannot be directly used for the thermalisation by simulating the corresponding Lindblad equation. Working out the Lindblad equation of the thermalisation requires the spectrum of the system BreuerPetruccione , which is the information that we want to obtain in the computation. To computation the thermal state, we have to assume that the Lindblad equation is unknown. Therefore, we need an environment to simulate the thermalisation. In this paper, we propose a hardware-efficient quantum algorithm for the simulation of Markovian and non-Markovian open-system dynamics, which can be used for solving the thermalisation and ground-state problems.

Qubits are valuable resources in the present and future. It is similar to the classical computational resources we use today but more severe. Fault-tolerant quantum computation based on the quantum error correction is the way to implement large-scale quantum computations, in which encoding one logical qubit may need thousands of physical qubits Fowler2012 ; OGorman2017 . Therefore, reducing the number of logical qubits is essential. Variational quantum algorithms for solving the ground-state problem or simulating the real and imaginary time evolutions have been developed recently Peruzzo2014 ; Wecker2015 ; Li2017 ; McArdle2018 , which can avoid the enormous qubit cost and are suitable for the near-term quantum computation. In this paper, our algorithm is in the category of conventional quantum algorithms demanding fault tolerance but does not rely on a good variational ansatz. We reduce the qubit cost by using a small environment to simulate the open-system dynamics induced by a big environment. We achieve it by reproducing reservoir correlation functions of the big environment in the small environment.

Open-system dynamics is determined by reservoir correlation functions. According to the expansion of the time-convolutionless (TCL) master equation, the simulation of open-system dynamics is more accurate if higher-order correlation functions are reproduced. By reproducing two-time correlation functions, we can simulate the Redfield equation and, therefore, the Markovian quantum master equation when the Markov approximation is justified. From the Markovian quantum master equation, we can simulate the thermalisation.

Our algorithm is beyond two-time correlation functions. Any $n$-time correlation functions can be reproduced, therefore we can simulate TCL master equation up to any $n$-th-order expansion. Reservoir correlation functions can be reproduced using tensor network, in which way the dimension of the environment increases exponentially with the number of terms in the system-environment coupling Luchnikov2019 . Our algorithm uses a different approach. By minimising the Hilbert space dimension for reproducing given correlation functions, the number of qubits required for representing the environment is $\sim\frac{n}{2}\log_{2}(N_{\omega}N_{\beta})$, where $N_{\omega}$ is the number of frequencies in the discretised environment spectrum, and $N_{\beta}$ is the number of terms in the coupling. Reservoir correlation functions can be exactly reproduced up to the spectrum discretisation, which usually converges polynomially with $N_{\omega}$.

The theory of open-system dynamics is introduced in Sec. II. To simulate the open system dynamics given by the Hamiltonian $H$ and the environment state $\rho_{\rm E}$, instead, we implement the dynamics of the Hamiltonian $\widetilde{H}$ and the environment state $\widetilde{\rho}_{\rm E}$ on the quantum computer, as shown in Fig. 1. The overview of the algorithm is given in Sec. III, and details of the algorithm are discussed in Sec. IV, V and VI. In Sec. VII, we discuss how to reinitialise the environment in the simulation. In Sec. VIII, the circuit implementation, qubit cost and gate-number cost are discussed. In Sec. IX, we give an illustrative example, and we numerically implement our algorithm on a classical computer to simulate the thermalisation of a qubit.

Later we will show how to choose $\widetilde{H}$ and $\widetilde{\rho}_{\rm E}$ such that reservoir correlation functions of $H$ and $\rho_{\rm E}$ can be reproduced on the quantum computer. Because we want to minimise the number of qubits representing the environment, the state of the environment may significantly change with time. Therefore, we may need to re-initialise the environment state during the simulation. We also show how to implement the re-initialise without significantly modifying correlation functions.

II Dynamics of open quantum systems weakly coupled to the environment

Given the Hamiltonian of the system and environment $H=H_{\rm S}+H_{\rm E}+\alpha H_{\rm I}$, where $H_{\rm S}$, $H_{\rm E}$ and $H_{\rm I}$ respectively denote Hamiltonians of the system, environment and interaction [see Fig. 1(a)], the evolution equation in the interaction picture is

[TABLE]

Here, $\rho$ is the state of the system and environment, $\alpha$ is a dimensionless coupling constant, and we have taken $\hbar=1$. Derived from this evolution equation, TCL equation BreuerPetruccione is

[TABLE]

for any initial state in the form $\rho(0)=\rho_{\rm S}(0)\otimes\rho_{\rm E}$, where $\mathcal{P}$ is a superoperator projection defined by $\mathcal{P}\rho\equiv\mathrm{Tr}_{\rm E}(\rho)\otimes\rho_{\rm E}$. We focus on the case that $\rho_{\rm E}$ is a stationary state of the environment, i.e. $[H_{\rm E},\rho_{\rm E}]=0$.

TCL master equation is the evolution equation of the system state $\rho_{\rm S}=\mathrm{Tr}_{\rm E}(\rho)$, because $\mathcal{P}\rho(t)=\rho_{\rm S}(t)\otimes\rho_{\rm E}$, in which only the system state $\rho_{\rm S}(t)$ evolves with time, and the environment state $\rho_{\rm E}$ is constant.

When the coupling between the system and environment is weak, the expansion of TCL generator $\mathcal{K}(t)$ in powers of the coupling constant $\alpha$ provides a series of approximate evolution equations. The expansion reads $\mathcal{K}(t)=\sum_{n=1}^{\infty}\alpha^{n}\mathcal{K}_{n}(t)$, where $\mathcal{K}_{n}(t)$ does not depend on $\alpha$. For example, up to the fourth order, we have $\mathcal{K}_{1}(t)=\mathcal{P}\mathcal{L}(t)\mathcal{P}$, $\mathcal{K}_{2}(t)=\int_{0}^{t}dt_{1}\mathcal{P}\mathcal{L}(t)\mathcal{L}(t_{1})\mathcal{P}$, $\mathcal{K}_{3}(t)=\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\mathcal{P}\mathcal{L}(t)\mathcal{L}(t_{1})\mathcal{L}(t_{2})\mathcal{P}$ and

[TABLE]

where

[TABLE]

Without loss of generality, we have assumed $\mathcal{P}\mathcal{L}(t)\mathcal{P}=0$ for simplification. In general, the $n$-th-order TCL generator $\mathcal{K}_{n}(t)$ is determined by superoperators $\boldsymbol{\kappa}_{m}$ with $m\leq n$.

TCL equation is the exact evolution equation of the system, therefore, describes the non-Markovian dynamics of the system. By neglecting high-order terms and taking the approximation $\mathcal{K}_{2}(t)\simeq\mathcal{K}_{2}(\infty)$ under the assumption that the correlation time is short, we can get the Markovian quantum master equation $\frac{\partial}{\partial t}\mathcal{P}\rho(t)=\alpha^{2}\mathcal{K}_{2}(\infty)\mathcal{P}\rho(t)$.

Correlation functions of the environment

Open-system dynamics is determined by reservoir correlation functions. In general, the $n$-th-order TCL generator $\mathcal{K}_{n}(t)$ is determined by up to the $n$-time correlations BreuerPetruccione . The interaction can always be expressed in the form $H_{\rm I}=\sum_{\beta}A_{\beta}\otimes B_{\beta}$, where $A_{\beta}$ acts on the system, $B_{\beta}$ acts on the environment, and they are both Hermitian. Expanding $\boldsymbol{\kappa}_{n}$ using the expression of $H_{\rm I}$, we have

[TABLE]

Here, $\nu,\ldots,\nu_{n-1}=0,1$ are binary numbers indicating on which side the Hamiltonian acts, we define superoperators $\mathcal{C}_{\beta}(\iota,\nu)\rho\equiv[C_{\beta}(\iota)]^{\nu}\rho[C_{\beta}(\iota)]^{1-\nu}$, $C=A,B$, and $\iota=t,\omega$ ($\iota=\omega$ will be used later). Therefore, the superoperator $\boldsymbol{\kappa}_{n}$ is determined by system operators $A_{\beta}$ and $n$-time correlation functions of the environment. We note that reservoir correlation functions are time-ordered, i.e. $t\geq\cdots\geq t_{n-1}$.

For the environment with a discretised spectrum, the environment Hamiltonian can be decomposed according to the spectrum as $H_{\rm E}=\sum_{\varepsilon}\varepsilon\Pi(\varepsilon)$, where $\Pi(\varepsilon)$ is the projection onto the eigenspace of the eigenenergy $\varepsilon$. We define operators $B_{\beta}(\omega)\equiv\sum_{\varepsilon^{\prime}-\varepsilon=\omega}\Pi(\varepsilon)B_{\beta}\Pi(\varepsilon^{\prime})$, then $B_{\beta}(t)=\sum_{\omega}e^{-i\omega t}B_{\beta}(\omega)$. We can find that $B_{\beta}^{\dagger}(\omega)=B_{\beta}(-\omega)$, because $B_{\beta}$ is Hermitian.

Without loss of generality, we assume $\mathrm{Tr}(B_{\beta}\rho_{\rm E})=0$, i.e. $\mathcal{P}\mathcal{L}(t)\mathcal{P}=0$. We note that $\rho_{\rm E}$ is a stationary state. If $\mathrm{Tr}(B_{\beta}\rho_{\rm E})$ is not zero, we can replace $B_{\beta}$ with $B_{\beta}-\mathrm{Tr}(B_{\beta}\rho_{\rm E})\openone$ and $H_{\rm S}$ with $H_{\rm S}+\sum_{\beta}\mathrm{Tr}(B_{\beta}\rho_{\rm E})A_{\beta}\otimes\openone$, so that the total Hamiltonian is not changed but the assumption is satisfied.

Two distinct environments result in the same dynamics of the system if their correlation functions are the same and they are coupled to the system by the same set of operators $A_{\beta}$ [see Fig. 1(a)].

Theorem 1**.**

Let $\widetilde{H}=H_{\rm S}+\widetilde{H}_{\rm E}+\alpha\widetilde{H}_{\rm I}$ be the Hamiltonian of the system and a different environment, and $\widetilde{H}_{I}=\sum_{\beta}A_{\beta}\otimes\widetilde{B}_{\beta}$. The sufficient condition for the same dynamics of the system up to the $n$-th-order, i.e. $\sum_{m=1}^{n}\alpha^{m}\mathcal{K}_{m}(t)=\mathcal{P}\sum_{m=1}^{n}\alpha^{m}\widetilde{\mathcal{K}}_{m}(t)$, is that

[TABLE]

holds for all $m\leq n$.

III Overview of the algorithm

According to Theorem 1, in order to simulate the open-system dynamics driven by $H$ up to the $n$-th order, we can implement the evolution driven by $\widetilde{H}$ on the quantum computer. The algorithm has two stages. At the first stage, we compute correlation functions of the environment determined by $H_{\rm E}$ and $\{B_{\beta}\}$ and design the environment on the quantum computer, i.e. choose $\widetilde{H}_{\rm E}$, $\{\widetilde{B}_{\beta}\}$ and the dissipation of the environment to reproduce the same correlation functions. The purpose of our algorithm is to simulate the dynamics of an open quantum system and study the system rather than the environment. We assume that correlation functions of the environment are computable in classical computation. If correlation functions of the environment are classically intractable, we may need the quantum computer to study the dynamics of the environment, which is beyond the scope of this work. On the quantum computer, we want to minimise the size of the environment, therefore dissipation of the environment may be required in order to relax the environment and suppress the finite-size effect. We will give the protocol for designing the environment on the quantum computer later. At the second stage, we use the quantum computer to realise the time evolution driven by $\widetilde{H}$ and the dissipation [see Fig. 1(b)]. Given the corresponding Lindblad equation in the explicit form, the evolution can be realised on the quantum computer using a quantum circuit Bacon2001 ; Kliesch2011 ; Sweke2015 ; Candia2015 ; Sweke2016 ; Childs2017 ; Chenu2017 .

IV Simulation of the second-order equation and Markovian master equation

In this section, we consider the quantum simulation of the master equation with the second-order approximation. If $\mathcal{P}\mathcal{L}(t)\mathcal{P}=0$ and higher-order contributions are neglected, the evolution equation of the system reads $\frac{\partial}{\partial t}\mathcal{P}\rho(t)=\alpha^{2}\mathcal{K}_{2}(t)\mathcal{P}\rho(t)$, which can also be expressed in the from

[TABLE]

This equation is determined by two-time correlation functions

[TABLE]

In order to simulate the time evolution driven by Eq. (7), we reproduce such correlation functions on the quantum computer.

When the time scale over which reservoir correlation functions decay is negligible compared to the time scale over which the system evolves significantly, the Markov approximation is justified. Then $\mathcal{K}_{2}(t)\simeq\mathcal{K}_{2}(\infty)$ and the upper limit $t$ of the integral in Eq. (7) can be replaced by $\infty$. Our algorithm can simulate the open system dynamics with a finite correlation time, i.e. the dynamics is non-Markovian, but we focus on the case that the correlation time is short although may not be negligible.

IV.1 Algorithm for the second-order simulation

The correlation function of the environment can be expressed in the form

[TABLE]

where \gamma_{\beta,\beta^{\prime}}(\omega)=\mathrm{Tr}\big{[}B_{\beta}(\omega)B_{\beta^{\prime}}^{\dagger}(\omega)\rho_{\rm E}\big{]}. To choose the interaction operators $\{\widetilde{B}_{\beta}\}$ on the quantum computer, we diagonalise matrices $\gamma(\omega)$ on a classical computer. Matrices $\gamma(\omega)$ are Hermitian and positive. After the diagonalisation, we obtain $\gamma(\omega)=U(\omega)\Lambda(\omega)U^{\dagger}(\omega)$, where $\Lambda(\omega)$ is the diagonalised matrix, and $U(\omega)$ is unitary. Interaction operators $\{\widetilde{B}_{\beta}\}$ on the quantum computer depend on coefficients $g_{\beta,l}(\omega)=U_{\beta,l}(\omega)\sqrt{\Lambda_{l,l}(\omega)}$.

On the quantum computer, we use the Hilbert space $\widetilde{\mathcal{H}}_{\rm E}=\widetilde{\mathcal{H}}_{\rm v}\oplus\bigoplus_{\omega}\widetilde{\mathcal{H}}_{\omega}$ to represent the environment. Here, $\widetilde{\mathcal{H}}_{\rm v}$ is one-dimensional and contains only one state $|{\rm v}\rangle$ representing the vacuum, $\widetilde{\mathcal{H}}_{\omega}$ is $d_{\omega}$-dimensional and corresponds to the transition frequency $\omega$, and $d_{\omega}=\mathrm{rank}(\gamma(\omega))$. The orthonormal basis of $\widetilde{\mathcal{H}}_{\omega}$ is $\{|{\omega,l}\rangle\}$, where $l$ corresponds to the $l$-th eigenvalue of $\gamma(\omega)$. The dimension of the environment $\widetilde{\mathcal{H}}_{\rm E}$ is $d_{\rm E}=1+\sum_{\omega}d_{\omega}$, therefore we can use $N_{\rm E}=\lceil\log_{2}d_{\rm E}\rceil$ qubits to simulate the environment. We have $d_{\omega}\leq N_{\beta}$ and $d_{\rm E}\leq 1+N_{\omega}N_{\beta}$, where $N_{\beta}\equiv|\{\beta\}|$ is the number of terms in the interaction Hamiltonian, $N_{\omega}\equiv|\{\omega\}|\leq N_{\varepsilon}^{2}$ is the number of transition frequencies, and $N_{\varepsilon}\equiv|\{\varepsilon\}|$ is the number of eigenenergies in the discretised spectrum of the environment.

To simulate the environment, we take $\widetilde{\rho}_{\rm E}=|\rm v\rangle\langle\rm v|$,

[TABLE]

where $\sigma_{l}(\omega)\equiv|{\rm v}\rangle\langle\omega,l|$. Then, correlations functions can be reproduced on the quantum computer. We note that \mathrm{Tr}\big{[}\widetilde{B}_{\beta}(\omega)\widetilde{B}_{\beta^{\prime}}^{\dagger}(\omega)\widetilde{\rho}_{\rm E}\big{]}=[U(\omega)\sqrt{\Lambda(\omega)}\sqrt{\Lambda(\omega)}U^{\dagger}(\omega)]_{\beta,\beta^{\prime}}=\gamma_{\beta,\beta^{\prime}}(\omega), where

[TABLE]

Therefore, $\langle\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(t-s)\rangle=\langle B_{\beta}(t)B_{\beta^{\prime}}(t-s)\rangle$, where \langle\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(t-s)\rangle=\mathrm{Tr}\big{[}\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(t-s)\widetilde{\rho}_{\rm E}\big{]}.

IV.2 Discussion

In the algorithm, the initial state of the environment on the quantum computer is always the pure state $|{\rm v}\rangle$, and the pure state is not the ground state, because the frequency $\omega$ can take negative values (see Fig. 2). The system can release energy into the environment via a transition from the state $|{\rm v}\rangle$ to the state $|{{\omega,l}}\rangle$ with a positive $\omega$. Similarly, the system can absorb energy from the environment via a transition from the state $|{\rm v}\rangle$ to the state $|{{\omega,l}}\rangle$ with a negative $\omega$. For a thermal bath with the temperature $T$, the $\gamma$ matrix satisfies $\gamma_{\beta,\beta^{\prime}}(-\omega)=\exp(-\hbar\omega/k_{\rm B}T)\gamma_{\beta^{\prime},\beta}(\omega)$ BreuerPetruccione . If the temperature is [math], $\gamma(\omega)=0$, i.e. $g_{\beta,l}(\omega)=0$, for all negative $\omega$. Then, states $|{{\omega,l}}\rangle$ with a negative $\omega$ are decoupled from the system. In this case, the state $|{\rm v}\rangle$ is the effective ground state. If the temperature is finite, the system is not only coupled to positive-$\omega$ states but also negative-$\omega$ states. In this way, we can simulate a finite temperature environment using a pure state as the initial state of the environment.

V State space relevant to the $n$-th-order expansion

We can generalise the algorithm in Sec. IV to simulate the TCL master equation up to the $n$-th order. Before giving the general algorithm, we first analyse the space of relevant environment states. The dimension of the state space has an upper bound $d_{n,{\rm max}}=[(N_{\omega}N_{\beta})^{\lfloor n/2\rfloor+1}-1]/[N_{\omega}N_{\beta}-1]$ as we explain in the next paragraph. Given the dimension of the state space, we can use a Hilbert space with the same dimension as the environment on the quantum computer to simulate the $n$-th-order TCL master equation.

Now we explain the upper bound of the dimension. We can rewrite the $m$-time correlation function as

[TABLE]

Let $|{\psi}\rangle$ be the purification of the state $\rho_{\rm E}$, i.e. $|{\psi}\rangle$ is a state on the Hilbert space $\mathcal{H}_{\rm E}\otimes\mathcal{H}_{\rm a}$ satisfying $\mathrm{Tr}_{\rm a}(|\psi\rangle\langle\psi|)=\rho_{\rm E}$. Here, $\mathcal{H}_{\rm E}$ is the Hilbert space of the environment, and $\mathcal{H}_{\rm a}$ is the Hilbert space of an ancillary system with the minimum dimension $\mathrm{rank}(\rho_{\rm E})$. We define $b_{\beta}(\omega)\equiv B_{\beta}(\omega)\otimes\openone_{\rm a}$, where $\openone_{\rm a}$ is the identity operator of $\mathcal{H}_{\rm a}$. Expressing superoperators $\mathcal{B}$ using operators $B$, the last line of Eq. (13) can be expressed in the form

[TABLE]

and we can rewrite it as

[TABLE]

Then using the cyclic property of trace, we rewrite correlation functions in the form

[TABLE]

To simulate the $n$-th-order TCL master equation, we only need to consider correlation functions with $m\leq n$. Therefore, the maximum number of $b$ operators in the above equation is $n$. We introduce states

[TABLE]

where $\Omega=\omega+\omega_{1}+\cdots+\omega_{m-1}$. Then, correlation functions can always be expressed in the form

[TABLE]

where two arguments of $\phi$, $\beta^{\prime}$, $\omega^{\prime}$ and $\nu^{\prime}$ on the second line depends on $\beta$, $\omega$ and $\nu$ one the first line. Here, we remark that $b_{\beta}^{\dagger}(\omega)=b_{\beta}(-\omega)$. Because the maximum number of $b$ operators is $n$, all correlations can be expressed in the above form with $m_{\rm L},m_{\rm R}\leq\lfloor n/2\rfloor$. The number of states $|{\phi_{\Omega,m}(\cdots)}\rangle$ is $(N_{\omega}N_{\beta})^{m}$, because there are $m$ operators $b$ acting on $|{\psi}\rangle$, and each operator $b$ has $N_{\omega}N_{\beta}$ options. Therefore, the total number of all relevant states $V_{n}=\{|{\phi_{\Omega,m}(\cdots)}\rangle~{}|~{}m\leq\lfloor n/2\rfloor\}$ is $1+N_{\omega}N_{\beta}+\cdots+(N_{\omega}N_{\beta})^{\lfloor n/2\rfloor}$, which is the upper bound of the space dimension.

We can decompose the space of relevant states according to the frequency. Because $\rho_{\rm E}$ is a stationary state, the correlation function in Eq. (16) is nonzero only if the summation of frequencies is zero, i.e. $\omega+\omega_{1}+\cdots+\omega_{m-1}=0$. Therefore, two states $|{\phi_{\Omega_{\rm R},m_{\rm R}}(\cdots)}\rangle$ and $|{\phi_{\Omega_{\rm L},m_{\rm L}}(\cdots)}\rangle$ are orthogonal if $\Omega_{\rm R}\neq\Omega_{\rm L}$. Then, the space of relevant states can be decomposed as $\mathcal{H}_{\rm E}^{\rm r}=\bigoplus_{\Omega}\mathcal{H}_{\Omega}$, where $\mathcal{H}_{\Omega}$ is the span of states $\{|{\phi_{\Omega,m}(\cdots)}\rangle\}$ with the frequency $\Omega$.

VI General algorithm for simulating the environment

The algorithm has two stages. At the first stage, we compute correlation functions of the environment and work out how to encode the environment on the quantum computer. At the second stage, we use the quantum computer to realise the time evolution driven by a Hamiltonian worked out at first stage.

VI.1 Classical computation

To simulate the dynamics of an open quantum system up to the $n$-th-order expansion of TCL equation, we compute correlation functions, $g_{\phi,\phi^{\prime}}=\langle\phi|\phi^{\prime}\rangle$ and $b_{\phi,\phi^{\prime}}=\langle{\phi}|b|{\phi^{\prime}}\rangle$, where $|{\phi}\rangle,|{\phi^{\prime}}\rangle\in V_{n}$ and $b\in\{b_{\beta}(\omega)\}$. These correlations functions are all in the form of the last line in Eq. (16).

Using the Gram matrix $g_{\phi,\phi^{\prime}}$ and Gram-Schmidt orthogonalisation (see Appendix), we can obtain a $d_{\rm E}$-dimensional representation of states $|{\phi}\rangle$ and operators $b$, where $d_{\rm E}=\mathrm{rank}(g)\leq d_{n,{\rm max}}$ is the dimension of $\mathcal{H}_{\rm E}^{\rm r}$. Each state $|{\phi}\rangle\in V_{n}$ maps to a $d_{\rm E}$-dimensional vector $|{\widetilde{\phi}}\rangle$, and each $b\in\{b_{\beta}(\omega)\}$ maps to a $d_{\rm E}$-dimensional matrix $\widetilde{b}$. These $d_{\rm E}$-dimensional vectors and matrices satisfy $\langle\widetilde{\phi}|\widetilde{\phi}^{\prime}\rangle=g_{\phi,\phi^{\prime}}$ and $\langle{\widetilde{\phi}}|\widetilde{b}|{\widetilde{\phi}^{\prime}}\rangle=b_{\phi,\phi^{\prime}}$. Then,

[TABLE]

holds for all $m$-th order correlation functions if $m\leq 2\lfloor n/2\rfloor+1$. Given $|{\widetilde{\phi}}\rangle$ and $\widetilde{b}$, we can simulate dynamics of the open quantum system on the quantum computer.

As the same as $\mathcal{H}_{\rm E}^{\rm r}$, the space of vectors $\{|{\widetilde{\phi}}\rangle\}$ can be decomposed in the form $\widetilde{\mathcal{H}}_{\rm E}=\bigoplus_{\Omega}\widetilde{\mathcal{H}}_{\Omega}$, where $\widetilde{\mathcal{H}}_{\Omega}$ is the span of states $\{|{\widetilde{\phi}_{\Omega,m}(\cdots)}\rangle\}$ with the frequency $\Omega$, because $\langle\widetilde{\phi}_{\Omega_{\rm L},m_{\rm L}}(\cdots)|\widetilde{\phi}_{\Omega_{\rm R},m_{\rm R}}(\cdots)\rangle=0$ if $\Omega_{\rm R}\neq\Omega_{\rm L}$. We remark that $|{\widetilde{\psi}}\rangle$ is in the subspace $\widetilde{\mathcal{H}}_{\Omega=0}$.

VI.2 Quantum computation

The simulation performed on the quantum computer is as follows. On the quantum computer, we use a $d_{\rm E}$-dimensional Hilbert space $\widetilde{\mathcal{H}}_{\rm E}=\bigoplus_{\Omega}\widetilde{\mathcal{H}}_{\Omega}$, i.e. $N_{\rm E}=\lceil\log_{2}d_{\rm E}\rceil$ qubits, to represent the environment, where $d_{\rm E}\leq d_{n,{\rm max}}$. We use $\widetilde{\Pi}_{\Omega}$ to denote the orthogonal projection on the subspace $\widetilde{\mathcal{H}}_{\Omega}$.

To simulate the environment, we take $\widetilde{\rho}_{\rm E}=|\widetilde{\psi}\rangle\langle\widetilde{\psi}|$,

[TABLE]

On the quantum computer, we implement the time evolution with the Hamiltonian $\widetilde{H}=H_{\rm S}+\widetilde{H}_{\rm E}+\alpha\widetilde{H}_{\rm I}$ and the environment initial state $\widetilde{\rho}_{\rm E}$. Then TCL generator of the system evolution on the quantum computer $\widetilde{\mathcal{K}}$ is the same as the generator of the dynamics to be simulated $\mathcal{K}$ up to the $n$-th-order expansion, i.e. $\widetilde{\mathcal{K}}_{m}(t)=\widetilde{\mathcal{P}}\mathcal{K}_{m}(t)$ for all $m\leq n$, according to Theorem 1. The proof is given in Appendix.

VI.3 Discussion

We can understand the algorithm as follows. By introducing the ancillary Hilbert space $\mathcal{H}_{\rm a}$, we can write the purification of the initial state $\rho_{\rm E}$ as $|{\psi}\rangle=\sum_{\varepsilon}\sqrt{p_{\varepsilon}}|{\Psi_{\varepsilon}}\rangle_{\rm E}\otimes|{\Phi_{\varepsilon}}\rangle_{\rm a}$, where $|{\Psi_{\varepsilon}}\rangle_{\rm E}$ is the eigenstate of the environment with the energy $\varepsilon$, and both $\{|{\Psi_{\varepsilon}}\rangle_{\rm E}\}$ and $\{|{\Phi_{\varepsilon}}\rangle_{\rm a}\}$ are orthonormal. Here, we have used that $\rho_{\rm E}$ is a stationary state. Then, we can write the Hamiltonian of the system, environment and ancillary system as $H^{\prime}=H_{\rm S}+H_{\rm E}^{\prime}+H_{\rm I}$, where $H_{\rm E}^{\prime}=H_{\rm E}+H_{\rm a}$ and $H_{\rm a}=-\sum_{\varepsilon}\varepsilon|\Phi_{\varepsilon}\rangle\langle\Phi_{\varepsilon}|$. According to $H^{\prime}$, the ancillary system is decoupled from the system and environment, and $|{\psi}\rangle$ is an eigenstate of $H_{\rm E}^{\prime}$ with the energy [math]. Let $\Pi_{\rm E}$ be the orthogonal projection onto the relevant subspace $\mathcal{H}_{\rm E}^{\rm r}$, then $\widetilde{H}=\Pi_{\rm E}H^{\prime}\Pi_{\rm E}$.

We would like to remark that, the ancillary system discussed here has been included in the environment $\widetilde{H}_{\rm E}$ on the quantum computer, which are not the ancillary qubits used for realising the evolution circuit shown in Fig. 1(b).

Similar to the second-order simulation, the initial state of the environment on the quantum computer is always a pure state, and the pure state is not the ground state, because the frequency $\Omega$ can take both positive and negative values (see Fig. 2).

VII Relaxation of the environment on the quantum computer

Usually, higher-order terms of TCL equation are less significant, because of not only the weak coupling but also the huge energy and information capacity of the environment, i.e. the influence of the system on the environment is small. However, on the quantum computer, the environment always has a finite size. As a result, high-order terms may become significant when the evolution time is long enough, specifically when the system and the environment exchange multiple excitations and the environment becomes saturate. Therefore, in this case we need to introduce the relaxation of the environment, i.e. the dynamics implemented on the quantum computer is modified to $\frac{\partial}{\partial t}\rho(t)=-i[\widetilde{H},\rho(t)]+\mathcal{L}_{\rm R}\rho(t)$, where the Lindblad superoperator $\mathcal{L}_{\rm R}$ acts on the environment and causes the relaxation. Evolution of such a Lindblad equation can also be implemented on the quantum computer Bacon2001 ; Kliesch2011 ; Sweke2015 ; Candia2015 ; Sweke2016 ; Childs2017 ; Chenu2017 . In this section, we present three protocols for the environment relaxation.

Before we give relaxation protocols, we take the algorithm for the Markovian master equation simulation in Sec. IV as an example to show the impact of the finite environment. According to $\widetilde{H}$, we have $\widetilde{\mathcal{K}}_{1}(t)=\widetilde{\mathcal{K}}_{3}(t)=0$, and $\widetilde{\mathcal{K}}_{4}(t)$ has four terms as shown in Eq. (3). The condition of the Markov approximation is the short correlation time $\tau_{\rm E}$ of the environment, i.e. $\langle B_{\beta}(t)B_{\beta^{\prime}}(t-s)\rangle$ is insignificant if $s>\tau_{\rm E}$. Then, $\boldsymbol{\kappa}_{2}(t,t-s)$ is insignificant if $s>\tau_{\rm E}$. As a result, integrals of the last two terms in $\mathcal{K}_{4}(t)$ leads to $\mathcal{O}(\tau_{\rm E}^{3})$. For example, the term $\boldsymbol{\kappa}_{2}(t,t_{2})\boldsymbol{\kappa}_{2}(t_{1},t_{3})$ is significant only in the region defined by $t\geq t_{1}\geq t_{2}$, $t\geq t_{2}\geq t-\tau_{\rm E}$ and $t_{2}\geq t_{3}\geq t_{1}-\tau_{\rm E}$. It is similar for $\boldsymbol{\kappa}_{2}(t,t_{3})\boldsymbol{\kappa}_{2}(t_{1},t_{2})$. However, integrals of the second term result in $\mathcal{O}(\tau_{\rm E}^{2}t)$, because $\boldsymbol{\kappa}_{2}(t,t_{1})\boldsymbol{\kappa}_{2}(t_{2},t_{3})$ is significant if $t-t_{1}\leq\tau_{\rm E}$, $t_{2}-t_{3}\leq\tau_{\rm E}$, but $t_{1}-t_{2}$ can be any value. Therefore, $\mathcal{K}_{4}(t)$ is small only if the second term and the first term $\boldsymbol{\kappa}_{4}(t,t_{1},t_{2},t_{3})$ cancel with each other, i.e. $\boldsymbol{\kappa}_{4}(t,t_{1},t_{2},t_{3})\simeq\boldsymbol{\kappa}_{2}(t,t_{1})\boldsymbol{\kappa}_{2}(t_{2},t_{3})$ when $t_{1}-t_{2}>\tau_{\rm E}$, which means that two excitations in the environment do not interfere with each other if they are separated by a time interval bigger than $\tau_{\rm E}$. However, in our algorithm for simulating the second-order equation, at most only one excitation can exist in the environment on the quantum computer, and the first excitation always prevents the second excitation, therefore $\widetilde{\mathcal{K}}_{4}(t)=\mathcal{O}(\tau_{\rm E}^{2}t)$.

As an example, we consider one of sixteen of terms in $\widetilde{\boldsymbol{\kappa}}_{4}(t,t_{1},t_{2},t_{3})$,

[TABLE]

Because at most only one excitation can exist in the environment, the contribution of the following components is nonzero [see Eq. (11)]: the $\sigma_{l}^{\dagger}(\omega)$ component of $\widetilde{B}_{\beta_{3}}$, the $\sigma_{l}(\omega)$ component of $\widetilde{B}_{\beta_{2}}$, the $\sigma_{l}^{\dagger}(\omega)$ component of $\widetilde{B}_{\beta_{1}}$ and the $\sigma_{l}(\omega)$ component of $\widetilde{B}_{\beta}$. As a result,

[TABLE]

which is significant if $t-t_{3}\leq\tau_{\rm E}$. We remark that $t_{1}$ and $t_{2}$ are between $t$ and $t_{3}$. The corresponding term in $\widetilde{\boldsymbol{\kappa}}_{2}(t,t_{1})\widetilde{\boldsymbol{\kappa}}_{2}(t_{2},t_{3})$ is

[TABLE]

where

[TABLE]

For any value of $t_{1}-t_{2}$, the corresponding term in $\widetilde{\boldsymbol{\kappa}}_{2}(t,t_{1})\widetilde{\boldsymbol{\kappa}}_{2}(t_{2},t_{3})$ can be significant. Therefore, $\widetilde{\boldsymbol{\kappa}}_{4}(t,t_{1},t_{2},t_{3})\simeq\widetilde{\boldsymbol{\kappa}}_{2}(t,t_{1})\widetilde{\boldsymbol{\kappa}}_{2}(t_{2},t_{3})$ does not hold when $t_{1}-t_{2}>\tau_{\rm E}$.

Next, we show that $\widetilde{\mathcal{K}}_{4}(t)$ can be suppressed by introducing the environment relaxation.

VII.1 Reinitialisation protocol

A way to realise the environment relaxation is the periodic reinitialisation of the environment state at time $j\tau$, where $\tau$ is the period, and $j$ is an integer Terhal2000 ; Wang2011 . In such a protocol, correlation functions on the quantum computer with $m\leq n$ are significantly modified by the relaxation and become

[TABLE]

where $\mathcal{V}_{i}=[\openone]$ if $t_{i-1}$ ($t_{0}=t$) and $t_{i}$ are in the same period, i.e. $(j+1)\tau>t_{i-1}\geq t_{i}>j\tau$ for any integer $j$, otherwise $\mathcal{V}_{i}=\mathcal{P}$. Here, $\mathcal{P}$ is the projection onto the state $\rho_{\rm E}$, $\openone$ is the identity operator, and $[U]\rho=U\rho U^{\dagger}$.

For two-time correlation functions, $\langle\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(s)\rangle=\langle B_{\beta}(t)B_{\beta^{\prime}}(s)\rangle$ only if $t$ and $s$ are in the same period, otherwise it is zero. We note that even if $t$ and $s$ are close, the two-time correlation function is zero if they are in different periods. Because of the reinitialisation, $\widetilde{\boldsymbol{\kappa}}_{4}(t,t_{1},t_{2},t_{3})=\widetilde{\boldsymbol{\kappa}}_{2}(t,t_{1})\widetilde{\boldsymbol{\kappa}}_{2}(t_{2},t_{3})$ if $t_{1}-t_{2}>\tau$, therefore the fourth order term is suppressed to $\widetilde{\mathcal{K}}_{4}(t)=\mathcal{O}(\tau_{\rm E}^{2}\tau)$.

VII.2 Projective dissipation protocol

We can implement the environment reinitialisation stochastically at a constant rate of $\Gamma$, i.e. the corresponding Lindblad superoperator is $\mathcal{L}_{\rm R}=\Gamma(\widetilde{\mathcal{P}}-[\openone])$. With such a relaxation term, correlation functions on the quantum computer with $m\leq n$ can also be expressed in the form of Eq. (23), but $\mathcal{V}_{i}=e^{-\Gamma s_{i}}[\openone]+(1-e^{-\Gamma s_{i}})\mathcal{P}$, where $s_{i}=t_{i-1}-t_{i}$. We remark that, $\mathcal{L}_{\rm R}$ and the environment Hamiltonian $-i[\widetilde{H}_{\rm E},\bullet]$ are commutative, because $\widetilde{\rho}_{E}$ is a stationary state.

Using the projective dissipation protocol, two-time correlation functions become $\langle\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(s)\rangle=e^{-\Gamma(t-s)}\langle B_{\beta}(t)B_{\beta^{\prime}}(s)\rangle$.

In some cases, correlation functions can be exactly reproduced on the quantum computer even with the presence of environment dissipation $\mathcal{L}_{\rm R}$. For the second-order equation simulation, if Fourier transformations of $e^{\Gamma(t-s)}\langle B_{\beta}(t)B_{\beta^{\prime}}(s)\rangle$ yield a set of positive matrices $\gamma(\omega)$, we can choose coefficients $g_{\beta,l}(\omega)$ according to $e^{\Gamma(t-s)}\langle B_{\beta}(t)B_{\beta^{\prime}}(s)\rangle$, so that $\langle\widetilde{B}_{\beta}(t)\widetilde{B}_{\beta^{\prime}}(s)\rangle=\langle B_{\beta}(t)B_{\beta^{\prime}}(s)\rangle$ when the dissipation is introduced. It is similar for higher-order equation simulations.

If correlation functions cannot be exactly reproduced, we may need to take $\Gamma\ll\tau_{\rm E}^{-1}$, so that correlation functions are not significantly modified. The relaxation time of the environment is $\Gamma^{-1}$. Therefore excitations in the environment do not interfere with each other if they are separated by a time interval bigger than $\Gamma^{-1}$, i.e. $\widetilde{\mathcal{K}}_{4}(t)=\mathcal{O}(\tau_{\rm E}^{2}\Gamma^{-1})$ in the second-order equation simulation.

VII.3 Conditional projective dissipation protocol

An optimal dissipation protocol relaxes the environment as soon as possible but does not modify correlations functions. Here we present such an environment dissipation protocol motivated by a typical open quantum system, an atom coupled to the free space in the vacuum state as shown in Fig. 3(a). The excited state of the atom decays into the ground state by emitting a photon into the free space. Because once the photon is emitted it leaves the atom and never comes back, the decay is irreversible. The correlation time depends on the length of the photon wave-package, because once the wave-package is out of the reach of the coupling, the photon cannot affect the atom anymore. Therefore, if a black material is placed at a finite but sufficient distance from the atom and absorbs photons [see Fig. 3(b)], the black material does not affect the evolution of the atom (neglecting the radiation from the material). We are interested in cases that the correlation time is short compared with the coupling between the system and environment, so that the expansion of TCL master equation is reasonable. We find that when the correlation time is short, only a subspace of the correlation-relevant state space $\mathcal{H}_{\rm E}^{\rm r}$ effectively contributes to correlation functions. Therefore, we can let the environment evolve without dissipation within the subspace, i.e. the left side of the black material, and the environment dissipates once its state is out of the subspace, i.e. the right side of the black material.

First, we consider the second-order equation simulation. We will generalise the protocol to higher-order equation simulations later. For the second-order equation simulation, we show that using the conditional dissipation protocol, the environment relaxes in the time scale $\sim\tau_{\rm E}$, but two-time correlation functions are only slightly modified.

Suppose that the environment spectrum is discretised with the uniform spacing $\delta\omega$, then each frequency $\omega$ corresponds to an integer $k$ and $\omega=k\delta\omega$. We apply the Fourier transformation to states $|{\omega,l}\rangle$ and define $|{x,l}\rangle\equiv\frac{1}{\sqrt{N_{\omega}}}\sum_{\omega}e^{-i\frac{\omega x}{c}}|{\omega,l}\rangle$, where $c=\frac{N_{\omega}\delta_{\omega}}{2\pi}$ and $x=0,\ldots,N_{\omega}-1$. These states form a ring as shown in Fig. 4. For a wave-package in the from $\sum_{x}a_{x}|{x,l}\rangle$, the evolution driven by $\widetilde{H}_{\rm E}$ transports the wave-package along the ring [see Fig. 4(a)], i.e. $e^{-i\widetilde{H}_{\rm E}t}(\sum_{x}a_{x}|{x,l}\rangle)=\sum_{x}a_{x}|{x+ct,l}\rangle$ when $ct$ is an integer. Here $|{x+N_{\omega}}\rangle\equiv|{x}\rangle$. Therefore, the evolution is periodic, and the period is $\frac{N_{\omega}}{c}=\frac{2\pi}{\delta\omega}$.

We would like to note that using the uniformly discretised spectrum on the quantum computer, two-time correlation functions with $s=t-t_{1}$ in the interval $[0,\frac{2\pi}{\delta\omega}]$ are reproduced in the form of Fourier series, which converges as $N_{\omega}\rightarrow\infty$. The optimal range of $k$ depends on correlation functions. Without loss of generality, we suppose $N_{\omega}$ is odd, and we take $k=-\frac{N_{\omega}-1}{2},\ldots,0,1,\ldots,\frac{N_{\omega}-1}{2}$.

In the $x$ representation, we can re-express interaction operators as

[TABLE]

where $\sigma_{l}(x)\equiv|v\rangle\langle x,l|$ and

[TABLE]

Therefore, $\widetilde{B}_{\beta}|{v}\rangle=\sum_{x,l}g_{\beta,l}^{*}(x)|{x,l}\rangle$ is a wave-package in the $x$ space. Without the dissipation, the correlation function

[TABLE]

is the overlap between two wavepackages $\widetilde{B}_{\beta}|{v}\rangle$ and $e^{-i\widetilde{H}_{\rm E}s}\widetilde{B}_{\beta^{\prime}}|{\rm v}\rangle$.

If $\beta=\beta^{\prime}$, the correlation function is maximised at $s=0$. The second wavepackage moves in the $x$ space with the speed $c$ without dispersion. As a result, the correlation function decreases with the time $s$. The correlation function vanishes at $s\sim\tau_{\rm E}$, which implies that the wavepackage is localised in the $x$ space with the width $~{}x_{\rm E}\equiv c\tau_{\rm E}$. Because the wavepackage is created by the coupling, the coupling strength $g_{\beta,l}(x)$ is also localised in the $x$ space with the same width as shown in Fig. 4. The localised coupling means that the matrix $\gamma(\omega)$ varies slowly with the frequency $\omega$. In the following, we assume that the coupling is localised in the region $0\leq x\leq x_{\rm E}$, which is reasonable when the system is coupled to the environment via local interactions.

The conditional dissipation protocol works as follows. In the region $0\leq x\leq x_{\rm T}$, where $x_{\rm T}\geq x_{\rm E}$, a wavepackage propagates freely without dissipation, such that two-time correlation functions can be reproduced. We remark that two-time correlation functions are only determined by the wavepackage in the region $0\leq x\leq x_{\rm E}$. In the region $x>x_{\rm T}$, the excitation decays at the rate of $\Gamma$, and the environment is stochastically reinitialised to the state $|{\rm v}\rangle$.

To implement the conditional dissipation, at the rate of $\Gamma$ we perform a measurement to find out whether the environment is in states with $x>x_{\rm T}$, i.e. the projection $\Pi=\sum_{x>x_{\rm T},l}\sigma_{l}(x)^{\dagger}\sigma_{l}(x)$. The environment reinitialisation is implemented depending on the measurement outcome. The corresponding Lindblad superoperator reads $\mathcal{L}_{\rm R}=\Gamma(\widetilde{\mathcal{P}}[\Pi]+[\openone-\Pi]-[\openone])$. Because of the dissipation, the wavepackage disappears before the revival. When the wavepackage disappears, the environment is reinitialised, and the next excitation can enter the environment.

The dissipation may cause quantum Zeno effect, which can prevent the wavepackage from entering the dissipation region $x>x_{\rm T}$. The propagation from $|{x}\rangle$ to $|{x+1}\rangle$ takes the time $c^{-1}=\frac{2\pi}{N_{\omega}\delta\omega}$. Therefore, the quantum Zeno effect is weak if $\Gamma\ll c$. Here, $c^{-1}$ corresponds to the time resolution of the environment. When the time resolution is fine, we should have $\tau_{\rm E}\gg c^{-1}$. In this case, we can take $\Gamma=\tau_{\rm E}^{-1}$ and $x_{\rm T}=x_{\rm E}$, such that excitations in the environment do not interfere with each other if they are separated by a time interval bigger than $\tau_{\rm E}$, i.e. $\mathcal{K}_{4}(t)=\mathcal{O}(\tau_{\rm E}^{3})$.

As an example, let us consider a simple case that the interaction Hamiltonian of the system and environment consists of only one term, reading $\widetilde{H}_{\rm I}=\alpha A\otimes\widetilde{B}$, with

[TABLE]

where the environment spectrum is discretised with the uniform spacing $\delta\omega$, and $g(\omega)=a\sqrt{\frac{2\bar{\gamma}}{\omega^{2}+{\bar{\gamma}}^{2}}\frac{\delta\omega}{2\pi}}$, where $a$ and $\bar{\gamma}$ are constants with the dimensions of frequency. With such an environment, the reconstructed correlation function is $\langle\widetilde{B}(t)\widetilde{B}(t-s)\rangle=a^{2}e^{-\bar{\gamma}|s|}$ in the limit $N_{\omega}\rightarrow\infty$. In Fig. 5, we plot the wavepackage $\widetilde{B}|{\rm v}\rangle$ in the $x$ space with the conditional dissipation. One can see that the wavepackage travels freely from left to right, until it enters the dissipation zone in which it quickly diminishes.

Conditional reinitialisation. To avoid the quantum Zeno effect, we can replace the continuous-time dissipation with periodic conditional reinitialisation. The Lindblad superoperator becomes time-dependent and reads $\mathcal{L}_{\rm R}=\Gamma(t)(\widetilde{\mathcal{P}}[\Pi]+[\openone-\Pi]-[\openone])$, where $\Gamma(t)=+\infty$ when $t=jc^{-1}$, $\Gamma(t)=0$ when $t\neq jc^{-1}$, and $j$ is an integer. In other words, a wavepackage can propagate freely in each time interval with the length $c^{-1}$, i.e. the wavepackage can propagate by one site in the $x$ space, and at each time $t=jc^{-1}$ the environment is conditionally reinitialised by implementing the operation $\widetilde{\mathcal{P}}[\Pi]+[\openone-\Pi]$. In this way, without affecting the propagation of the wavepackage in the region $x\leq x_{\rm E}$, the environment relaxes in the time scale $\tau_{\rm E}$. Here we take $x_{\rm T}=x_{\rm E}$.

Generalisation to the higher-order simulations

In this section, we discuss how to generalise the conditional dissipation (or reinitialisation) protocol to higher-order simulations. For environments similar to the case in Fig. 3, the system only affects the part of the environment close to it, its influence (excitations) propagates in the environment, and the part close to the system relaxes in the time scale $\tau_{\rm E}$, i.e. correlations

[TABLE]

are negligible when $t_{i}-t_{i+1}>\tau_{\rm E}$. We remark that because correlation functions are reproduced, it is the same for the environment on the quantum computer. For such environments, the reinitialisation operation $\mathcal{P}$ on $[e^{-iH_{\rm E}s}]\mathcal{B}_{\beta_{i+1}}(t_{i+1},\nu_{i+1})\cdots\rho_{\rm E}$ does not affect correlation functions. Here $s>\tau_{\rm E}$, and $[e^{-iH_{\rm E}s}]$ is the superoperator denoting the free evolution of the environment. To implement the condition dissipation, we need to find the proper projection $\Pi$ representing the space of states that the influence of the system has left the interaction region.

As an example, we consider the fourth-order simulation using the environment on the quantum computer given by

[TABLE]

and

[TABLE]

Here, $|{\rm v}\rangle$ denotes the vacuum state and the initial state of the environment, i.e. $\widetilde{\rho}_{\rm E}=|\rm v\rangle\langle\rm v|$, $|{\omega}\rangle$ denotes the state of one excitation with the frequency $\omega$, and $|{\omega_{1},\omega_{2}}\rangle$ denotes the state of two excitations with frequencies $\omega_{1}$ and $\omega_{2}$, respectively. By choosing coupling coefficients $g_{\beta,\omega}$ and $g_{\beta,\omega,\omega_{1},\omega_{2}}$, we can reproduce some reservoir correlation functions (see Appendix). The general algorithm for higher-order simulations is given in Sec. VI.

Correlation functions reproduced in the environment given by Eq. (29) and Eq. (30) are

[TABLE]

[TABLE]

and

[TABLE]

Similar to the second-order simulation, we suppose that the environment is discretised with the uniform spacing $\delta$, i.e. $\omega,\omega_{1},\omega_{2}=k\delta\omega$ and $k=-\frac{N_{\omega}-1}{2},\ldots,0,1,\ldots,\frac{N_{\omega}-1}{2}$. Then, by applying the Fourier transformation, we have

[TABLE]

In the $x$ representation,

[TABLE]

where

[TABLE]

If we only consider correlation functions at discretised times, i.e. $cs,cs_{1},cs_{2},cs_{3}$ are integers, two-time and four-time correlation functions are

[TABLE]

and

[TABLE]

Three-time correlation functions are zero.

Correlations in Eq. (28) are negligible if the system is only coupled to environment states $|{x}\rangle$ and $|{x_{1},x_{2}}\rangle$ with $x,x_{2}=0,\ldots,x_{\rm E}$, where $x_{\rm E}=c\tau_{\rm E}$. We have $\Delta=0$ if $g_{\beta,x}=0$ when $x>x_{\rm E}$, $g_{\beta,x,x_{1},x_{2}}=0$ when $x_{2}>x_{\rm E}$, and $g_{\beta,x,x_{1},x_{2}}=\delta_{x,x_{1}}g_{\beta,x_{2}}$ when $x,x_{1}>x_{\rm E}$. It is obvious for two-time correlation functions. For four-time correlation functions, considering values of $\nu$ [see Eq. (5)], there are $16$ of them, but only four of them are independent, which are

[TABLE]

where $t\geq t_{1}\geq t_{2}\geq t_{3}$. We can check that $\Delta=0$ for all of them. Here, we have assumed that $N_{\omega}\gg x_{\rm E}$. Therefore, to implement the conditional dissipation, we can take the projection $\Pi=\sum_{x>x_{\rm E}}|x\rangle\langle x|+\sum_{x_{1}}\sum_{x_{2}>x_{\rm E}}|x_{1},x_{2}\rangle\langle x_{1},x_{2}|$.

VIII Circuit implementation, time cost and hardware resource requirement

Given the Hamiltonian $\widetilde{H}$, an initial state of the system $\rho_{\rm S}(0)$, and the initial state of the environment $\widetilde{\rho}_{\rm E}$, we can implement the unitary dynamics $\rho(t)=e^{-i\widetilde{H}t}\rho_{\rm S}(0)\otimes\widetilde{\rho}_{\rm E}e^{i\widetilde{H}t}$ on the quantum computer. Then, $\rho(t)$ is a solution of the evolution equation $\frac{\partial}{\partial t}\widetilde{\mathcal{P}}\rho(t)=\widetilde{\mathcal{K}}(t)\widetilde{\mathcal{P}}\rho(t)$. According to discussions in the previous section, $\widetilde{\mathcal{K}}(t)$ and $\mathcal{K}(t)$ are the same up to the $n$-th-order expansion.

We can implement the dynamics of $\widetilde{H}$ using the Trotterisation algoirthm Lloyd1996 . Let $N_{\rm S}$ be the number of qubits representing the system, then the total number of qubits used in the simulation is $N_{\rm S}+N_{\rm E}$. System operators can always be expanded using Pauli operators as $H_{\rm S}=\sum_{\sigma\in S_{H}}f_{H,\sigma}\sigma$ and $A_{\beta}=\sum_{\sigma\in S_{\beta}}f_{\beta,\sigma}\sigma$. Here, $S_{H}$ and $S_{\beta}$ are subsets of $N_{\rm S}$-qubit Pauli operators. Similarly, environment operators can also be expanded using Pauli operators as $\widetilde{H}_{\rm E}=\sum_{\sigma\in E_{H}}h_{H,\sigma}\sigma$ and $\widetilde{B}_{\beta}=\sum_{\sigma\in E_{\beta}}h_{\beta,\sigma}\sigma$. Here, $E_{H}$ and $E_{\beta}$ are subsets of $N_{\rm E}$-qubit Pauli operators. Expansion coefficients $f$ and $h$ are all real, because these expanded operators are all Hermitian. Using Trotterisation, the evolution implemented on the quantum computer is

[TABLE]

where $N_{\rm T}$ is the number of Trotter steps, and each exponential of $k$-qubit Pauli operator can be implemented on the quantum computer with up to $2(k-1)$ CNOT gates and $2k+1$ single-qubit gates Whitfield2011 . Therefore, the total number of gates $N_{\rm G}$ is less than $[(4N_{\rm S}-1)|S_{H}|+(4N_{\rm E}-1)|E_{H}|+(4N_{\rm S}+4N_{\rm E}-1)\sum_{\beta}|S_{\beta}||E_{\beta}|]N_{\rm T}$.

The Trotter-Suzuki decomposition is approximate, and the difference between $U_{N_{\rm T}}$ and $e^{-i\widetilde{H}}$ is

[TABLE]

where $\|\bullet\|$ denotes the operator norm, and $N_{\rm terms}=|S_{H}|+|E_{H}|+\sum_{\beta}|S_{\beta}||E_{\beta}|$ is the number of terms in the Hamiltonian. We can prove that $\|\widetilde{B}_{\beta}\|\leq\|B_{\beta}\|$ (see Appendix), therefore the norm of the Hamiltonian has the upper bound

[TABLE]

Here, $\max\{|\omega|\}\sim\|H_{\rm E}\|$. However, usually it is sufficient to truncate the frequency at $\max\{|\omega|\}\sim\|H_{\rm S}\|$ when the coupling is weak.

Usually, for a Hamiltonian with local interactions, the number of terms in the Hamiltonian, i.e. each of $|S_{H}|$, $|S_{\beta}|$ and $N_{\beta}$, is a polynomial with respect to the system size $N_{\rm S}$.

The number of qubits required for simulating the environment is $N_{\rm E}\sim\lfloor\frac{n}{2}\rfloor\log_{2}(N_{\omega}N_{\beta})$, because $d_{\rm E}\sim(N_{\omega}N_{\beta})^{\lfloor\frac{n}{2}\rfloor}$. According to the maximum number of environment Pauli operators, we have $|E_{H}|,|E_{\beta}|\sim 4^{N_{\rm E}}\sim(N_{\omega}N_{\beta})^{n}$.

To implement the conditional dissipation, we may need to introduce only one more qubit for the measurement of $\Pi$, i.e. we can use the state $|{1}\rangle$ of the qubit to indicate the subspace. Because the conditional dissipation operation is performed at a low rate, the cost of gate number is small compared with the unitary evolution.

In summary, the simulation requires $N_{\rm E}\sim\lfloor\frac{n}{2}\rfloor\log_{2}(N_{\omega}N_{\beta})$ qubits to simulate the environment. The number of terms in the Hamiltonian is $N_{\rm terms}=\mathcal{O}(N_{\omega}^{n}N_{\beta}^{n})$. Then, we need the number of Trotter steps to be $N_{\rm T}\sim N_{\rm terms}\|\widetilde{H}\|^{2}t^{2}/\epsilon_{\rm Trotter}=\mathcal{O}(N_{\omega}^{n}N_{\beta}^{n+1})$. Therefore the total number of gates is $N_{\rm G}=\mathcal{O}(N_{\rm E}N_{\omega}^{2n}N_{\beta}^{2n+2})$. We note that a variety of methods have been devoloped to reduce the gate number in the Trotterisation algorithm Berry2007 ; Wiebe2010 ; Berry2015 , which could be applied in our case.

In our algorithm, the system size can easily exceed the environment size. For example, to simulate the quantum master equation with $n=2$, considering an environment with a million discretised frequencies ($N_{\omega}=10^{6}$) and a thousand interaction terms ($N_{\beta}=10^{3}$), we only need about $30$ qubits for encoding the environment, which is even smaller than the system in a non-trivial quantum simulation problem (with above $50$ qubits).

IX Thermalisation of a qubit on quantum computer

Let us consider the thermalisation of a qubit at the zero temperature and finite temperature. The system Hamiltonian is $\widetilde{H}_{\rm S}=-\frac{\Delta}{2}\sigma^{\rm z}$. The system is coupled to the environment via only one term, i.e. $\widetilde{H}_{\rm I}=\alpha A\otimes\widetilde{B}$, where $A=\sigma^{\rm x}$, and $\widetilde{B}$ is the same as in Eq. (27), but coupling coefficients are Ritschel2014

[TABLE]

Here, $\beta=1/k_{\rm B}T$ is the temperature, and $a$ and $\bar{\gamma}$ are constants with the dimensions of frequency. We take $\rho_{\rm S}(0)=|+\rangle\langle+|$ as the initial state of the qubit, and $|{+}\rangle=\frac{1}{\sqrt{2}}(|{0}\rangle+|{1}\rangle)$.

In Fig. 6, we plot the probability in the ground state $p_{\rm g}=\langle{0}|\rho_{\rm S}(t)|{0}\rangle$, where $\rho_{\rm S}(t)$ is the state of the qubit at time $t$. When the simulated environment is at the zero temperature, i.e. $\beta\rightarrow+\infty$, the qubit evolves into the ground state $|{0}\rangle$, i.e. $p_{\rm g}$ goes to $1$. When the temperature is finite, the probability approaches a finite value and coincides with the thermal distribution. For comparison, we also plotted the probability in the evolution driven by the corresponding Lindblad equation of the thermalisation BreuerPetruccione . The difference between the environment-simulation result and the Lindblad equation result is due to the discretisation of the environment spectrum and approximations used to derive the Lindblad equation, including neglecting high-order terms in TCL equation and the Markovian approximation.

X Conclusions

In this paper, we propose a hardware-efficient quantum algorithm to simulate the TCL master equation up to any finite order. It is achieved by reproducing reservoir correlation functions using a minimised Hilbert space. The number of qubits representing the environment is $\sim\lfloor\frac{n}{2}\rfloor\log_{2}(N_{\omega}N_{\beta})$ in the $n$-th-order simulation. We remark that $n=2$ in the simulation of the Markovian quantum master equation and the thermalisation. In our algorithm, the system size can easily exceed the environment size, e.g. when the system has tens of qubits. Because the environment on the quantum computer is small, it needs to be reinitialised in the simulation of a long-time evolution. We also propose an efficient reinitialisation protocol without significantly changing reservoir correlation functions. We illustrate our algorithm by using a classical computer and numerically simulate the thermalisation of a qubit at the zero and finite temperatures. Our results pave the way for practical quantum open-system simulation on a universal quantum computer.

Acknowledgements.

This work is supported by National Natural Science Foundation of China (Grant No. 11875050) and NSAF (Grant No. U1930403). HYS is also supported by China Postdoctoral Science Foundation (Grant No. 2018M630063) and National Natural Science Foundation of China (Grant No. 11905209).

Appendix A Gram-Schmidt orthogonalisation

In this section, we explicitly present the Gram-Schmidt orthogonalisation process. We have $d_{n,{\rm max}}$ vectors $V_{n}=\{|{\phi_{\Omega,m}(\cdots)}\rangle\}$ in $\mathcal{H}_{\rm E}^{\rm r}$, where $m\leq\lfloor n/2\rfloor$. We label these vectors as $|{\phi_{1}}\rangle,\ldots,|{\phi_{d_{n,{\rm max}}}}\rangle$. Without loss of generality, we take $|{\phi_{1}}\rangle=|{\psi}\rangle$, which can simplify the preparation of the environment initial state on the quantum computer, and we assume that states from $|{\phi_{1}}\rangle$ to $|{\phi_{d_{\rm E}}}\rangle$ are linearly independent. We note that $g$ is a $d_{n,{\rm max}}$-dimensional matrix with rank $d_{\rm E}$. The state $|{\psi}\rangle$ is normalised, therefore we take $|{e_{1}}\rangle=|{\phi_{1}}\rangle$. Then, we can obtain $d_{\rm E}$ orthonormal basis states by iterating

[TABLE]

Given $|{e_{j}}\rangle=\sum_{\phi\in V_{n}}e_{j,\phi}|{\phi}\rangle$, we compute the overlap using $\langle e_{j}|\phi_{i}\rangle=\sum_{\phi\in V_{n}}e_{j,\phi}^{*}g_{\phi,\phi_{i}}$. The outcome of the Gram-Schmidt orthogonalization is the $d_{\rm E}\times d_{n,{\rm max}}$ matrix $e_{j,\phi}$.

Using the matrix $e_{j,\phi}$, we can express states $|{\phi}\rangle\in V_{n}$ and operators $b\in\{b_{\beta}(\omega)\}$ using the orthonormal basis of the subspace $\mathcal{H}_{\rm E}^{\rm r}=\mathrm{span}(V_{n})$, i.e. $|{\phi}\rangle=\sum_{i}\langle e_{i}|\phi\rangle|{e_{i}}\rangle$ and $\Pi_{\rm E}b\Pi_{\rm E}=\sum_{i,j}\langle{e_{i}}|b|{e_{j}}\rangle|e_{i}\rangle\langle e_{j}|$, where $\Pi_{\rm E}=\sum_{i=1}^{d_{\rm E}}|e_{i}\rangle\langle e_{i}|$ is the projection onto the subspace, $\langle e_{i}|\phi\rangle=\sum_{\phi^{\prime}\in V_{n}}e_{i,\phi^{\prime}}^{*}g_{\phi^{\prime},\phi}$ and $\langle{e_{i}}|b|{e_{j}}\rangle=\sum_{\phi^{\prime},\phi\in V_{n}}e_{i,\phi^{\prime}}^{*}g_{\phi^{\prime},\phi}e_{j,\phi}$.

Let $\{|{\widetilde{e_{i}}}\rangle~{}|~{}i=1,\ldots,d_{\rm E}\}$ be $d_{\rm E}$-dimensional orthonormal states. Each $|{\widetilde{e_{i}}}\rangle$ is a state in $\widetilde{\mathcal{H}}_{\rm E}$ on the quantum computer. Then, for operators $b\in\{b_{\beta}(\omega)\}$, we define $\widetilde{b}\equiv\sum_{i,j}\langle{e_{i}}|b|{e_{j}}\rangle|\widetilde{e}_{i}\rangle\langle\widetilde{e}_{j}|$. Because $b_{\beta}^{\dagger}(\omega)=b_{\beta}(-\omega)$, $[\Pi_{\rm E}b_{\beta}(\omega)\Pi_{\rm E}]^{\dagger}=\Pi_{\rm E}b_{\beta}(-\omega)\Pi_{\rm E}$. Therefore, $\widetilde{b}_{\beta}^{\dagger}(\omega)=\widetilde{b}_{\beta}(-\omega)$.

For a state $|{\varphi}\rangle\in\mathrm{span}(V_{n})$, we define $|{\widetilde{\varphi}}\rangle\equiv\sum_{i}\langle e_{i}|\varphi\rangle|{\widetilde{e}_{i}}\rangle$. Then,

[TABLE]

We remark that $\Pi_{\rm E}|{\varphi}\rangle=|{\varphi}\rangle$.

Because $\Pi_{\rm E}b_{m}\Pi_{\rm E}\cdots\Pi_{\rm E}b_{2}\Pi_{\rm E}b_{1}|{\psi}\rangle=b_{m}\cdots b_{2}b_{1}|{\psi}\rangle\in V_{n}$ for all $m\leq\lfloor n/2\rfloor$, the following equation holds for all $m\leq 2\lfloor n/2\rfloor+1$,

[TABLE]

Therefore, Eq. (17) holds for all $m\leq 2\lfloor n/2\rfloor+1$.

Appendix B Proof of the algorithm

Because $\langle{\widetilde{\phi}_{\Omega_{\rm L},m_{\rm L}}(\cdots)}|\widetilde{b}_{\beta}(\omega)|{\widetilde{\phi}_{\Omega_{\rm R},m_{\rm R}}(\cdots)}\rangle=0$ if $\omega\neq\Omega_{\rm L}-\Omega_{\rm R}$, we have $\widetilde{b}_{\beta}(\omega)=\sum_{\Omega_{\rm L}-\Omega_{\rm R}=\omega}\widetilde{\Pi}_{\Omega_{\rm L}}\widetilde{b}_{\beta}(\omega)\widetilde{\Pi}_{\Omega_{\rm R}}$, then $\widetilde{B}_{\beta}(t)=\sum_{\omega}e^{-i\omega t}\widetilde{b}_{\beta}(\omega)$. Therefore, $\widetilde{B}_{\beta}(\omega)=\widetilde{b}_{\beta}(\omega)$.

Correlation functions on the quantum computer can be expressed as

[TABLE]

Because of Eq. (17), the following equation holds for all $m\leq n$,

[TABLE]

Therefore, Eq. (6) holds for all $m\leq n$.

Appendix C Norm of $\widetilde{B}_{\beta}$

For any state $|{\widetilde{\varphi}}\rangle\in\widetilde{\mathcal{H}}_{\rm E}$, we have

[TABLE]

where $b_{\beta}=\sum_{\omega}b_{\beta}(\omega)=B_{\beta}\otimes\openone_{\rm a}$. Notice that $\||{\varphi}\rangle\|=\||{\widetilde{\varphi}}\rangle\|$ and $\|b_{\beta}\|=\|B_{\beta}\|$, we have $\|\widetilde{B}_{\beta}\|\leq\|B_{\beta}\|$.

Appendix D Correlation functions reproduced in the fourth-order example

Correlation functions reproduced in the environment given by Eq. (29) and Eq. (30) are

[TABLE]

[TABLE]

and

[TABLE]