用于模拟非马尔可夫开放量子系统的辅助量子比特链量子算法

《Physical Review A》：Ancilla-train quantum algorithm for simulating non-Markovian open quantum systems

　　

我们提出了一种量子算法，用于模拟与高斯环境耦合的开放量子系统，该算法适用于任意构型和任意耦合强度。该算法可应用于强耦合或非马尔可夫环境的问题、多个环境之间未达热平衡的问题以及含时依赖哈密顿量的问题。我们证明，该算法能够以任意精度复现此类问题的真实动力学，并且对于广泛的问题范围，相对于Trotter化时间演化仅增加微小的资源开销；该开销关于目标精度的倒数呈低阶多项式。该算法基于一个核心洞察：任何高斯环境都可以表示为一列辅助量子比特，它们通过与系统的时间局域耦合依次相互作用，该耦合由浴关联函数的卷积平方根给出；这是我们的工作的一个次要结果。在自旋-玻色子模型的数值模拟中，我们发现，在中等至强耦合区域，使用适中数量（101–103）的辅助量子比特即可实现精确模拟。我们的成果为量子计算机在高效模拟非平衡和开放量子系统方面的应用打开了大门。

I. 引言

开发量子计算机的一个关键动力，在于其能够极大地加速由量子力学支配的物理过程的模拟[1]。为了发挥这一潜在能力，人们开发了多种量子算法来模拟自然界中的量子系统[2–11]。大量努力集中在孤立和平衡量子系统上，这些系统可以通过基态、激发谱或热系综来表征[5,9,12–22]。与此同时，自然界中许多重要的过程涉及开放和非平衡量子系统，它们具有动力学特征、结构化环境、时依赖力或处于相互非平衡状态的多个浴。非平衡和开放量子系统出现在生物学[23]、化学[24–26]或光捕获[27,28]等领域。虽然已有各种经典方法被开发用于解决这些问题[29–43]，但在经典计算机上模拟相对简单的系统仍然是具有挑战性的，例如由于非马尔可夫效应或热力学原理的失效[35,36,44]。这就引发了一个问题：我们如何有效地利用量子计算机来模拟开放和非平衡量子系统？

在本工作中，我们通过引入一种量子算法来回答上述问题，该算法能够高效地模拟与高斯环境耦合的任何开放量子系统的动力学和稳态。该算法被称为辅助量子比特列车算法（ATA），例如可应用于强耦合或非马尔可夫环境的问题、多个环境之间未达热平衡的问题以及时依赖哈密顿量的问题。ATA由标准的Trotter化演化构成，其中每次迭代都通过一个门序列将系统观测量与一个有限的辅助量子比特阵列——即一列辅助量子比特——耦合起来；见图1。通过适当但简单的系统-辅助列车耦合，我们证明ATA能够以任意精度高效模拟具有高斯环境的任何开放量子系统，精确复现由底层微观模型产生的真实动力学。

图1. 辅助量子比特列车算法（ATA）：ATA通过Trotter化演化，以任意精度模拟具有高斯环境的任何开放量子系统的演化。每次迭代——此处为δt——通过两个门序列将系统传播一个时间步δt，其中$U_S$是由系统哈密顿量产生的时间演化，而$U_{SA}$将一个或多个系统观测量与一个辅助量子比特寄存器（红色、紫色和蓝色框；黑色线条表示耦合）连接起来。迭代结束后，一部分辅助量子比特被带到下一个[$\delta t$]迭代（紫色），其余的被丢弃（红色），并在$|0\rangle$态中引入相同数量的新辅助量子比特（蓝色）。通过从浴的跃迁关联子[39]获得的适当系统-辅助耦合结构，我们证明ATA能够以有限辅助寄存器将系统的真实动力学复现至任意精度。

我们为ATA提供了严格的误差界和资源估计。特别是，我们证明对于广泛的一类问题，ATA相对于Trotter化哈密顿量模拟仅产生微小的开销：即严格局域或准局域的系统-辅助耦合电路（取决于模型类别），加上一个有限的辅助寄存器。作为另一个优势，我们认为ATA能够容忍有限的硬件错误率，因为其输出——开放量子系统的不动点轨迹或稳态——对通用噪声具有微扰鲁棒性。

我们基于一个核心洞察构建ATA：任何高斯环境都可以表示为一列量子比特——或量子二进制白噪声信号——通过适当的时间局域相互作用依次与系统交互。具体而言，我们证明当时间局域相互作用与浴关联函数（或跃迁关联子，在文献[39]中识别）的卷积平方根成正比时，这种表示是精确的。我们将这种表示称为高斯量子噪声的辅助量子比特列车表示（ATR）。ATR提供了一种系统的时域碰撞模型，以量子比特表示任何高斯量子环境，与频域赝模式玻色子表示[31,32,35,42]互补。作为关键的次要结果，我们获得了截断ATR辅助列车所产生误差的严格界，从而能够以有限量子比特寄存器系统地近似表示任何高斯环境。利用ATR所提供的高斯环境自然量子比特表示，我们通过对所得时间演化进行Trotter化，简单地获得了ATA。

有趣的是，近期在开放量子系统文献中出现了一个与ATR相关的构造[43]，它涉及一列谐振子而非量子比特，以及一个基于正交多项式的构造，不具有ATR的离散时间平移对称性。

ATR的一个说明性特例出现在对称功率谱密度情况下，此时环境等效于作用于系统的经典噪声信号[44]。在此情况下，ATR意味着任何具有高斯关联的经典噪声信号——即彩色噪声——都可以通过足够密集的随机符号相同脉冲序列生成；见图3。因此，ATR可以被视为经典结果——彩色噪声信号可以通过过滤（二进制）白噪声信号生成[45]——的量子噪声推广；此处跃迁关联子提供了滤波核。

ATA与近期为开放量子系统的系统量子模拟出现的另外两种方法互补：量子Gibbs采样（QGS）算法[12,17–22]和基于系统地推导出的非久塞林德方程的马尔可夫碰撞模型算法[39–41,46]，此处称为系统林德（SL）算法；详见第V A节（另见文献[47]，其中提出了一个先于这些发展的方案）。我们在图2和表I中总结了ATA、QGS和SL算法的适用范围和优势。与QGS和SL算法相比（后者通常仅在弱系统-环境耦合极限下才精确[48,49]，且需要编译准局域跃迁算符和Lamb位移），ATA提供以下关键优势：(1) ATA对系统-环境耦合的所有阶都精确，因此可以有效地捕获非马尔可夫效应。(2) ATA不需要编译准局域跃迁算符或Lamb位移（即系统哈密顿量的重整化），而QGS或SL算法需要此步骤，因此可能只需在每次Trotter迭代中增加一个低深度酉电路即可实现。(3) 与SL类似但不同于QGS，ATA可应用于非平衡问题，并能给出动力学、不动点轨迹以及稳态。ATA相对于QGS和SL算法的代价是一个有限的辅助量子比特寄存器，对于给定的固定目标精度，其与系统-浴耦合呈线性缩放。因此，当电路深度是限制资源时，ATA提供了广泛的适用性和高效实现。

表I. 系统开放量子系统模拟量子算法的比较。我们比较了ATA、量子Gibbs采样（QGS）和系统林德（SL）算法的特征；这些算法的定义见正文。此处$\epsilon$表示算法的目标精度，而$g$、$\tau_c$和$\Lambda$分别表示系统-浴耦合、关联时间和环境紫外截断的特征尺度，定义于式(17)和(20)中。

| 特性 | QGS | SL | ATA |
|---|---|---|---|
| 能否用准局域系统-辅助耦合实现 | ? | ? | ? |
| 严格局域系统-辅助耦合 | ? | ? | ? |
| 无需Lamb位移编译 | ? | ? | ? |
| 能否返回热态 | ? | (?) | (?) |
| 非平衡稳态 | ? | ? | ? |
| 动力学 | ? | ? | ? |
| 能否模拟非平衡浴 | ? | ? | ? |
| 时依赖哈密顿量 | ? | ? | ? |
| 非马尔可夫区域 | ? | ? | ? |
| 非高斯环境 | ? | ? | ? |
| 所需辅助量子比特数 | 1 | 1 | N |

a ATA在每次Trotter迭代中涉及严格局域辅助-系统耦合，条件是模拟中与环境耦合的观测量可以表示为系统寄存器中的局域观测量。

b 文献[19–21]发现了具有严格局域系统-辅助耦合和多个辅助量子比特的近似QGS类算法，这些算法可能产生近似的（非精确的）Gibbs态。

c 此处Lamb位移指哈密顿量的准局域重整化，对于QGS和SL算法对应的林德过程是必需的。

d 物理系统仅在系统-浴耦合趋零的极限下才收敛到Gibbs态[48]；(?)表示稳态包含对此事实时修正，而?表示算法返回精确的Gibbs态。

e 例如，对于基于文献[39,41]中林德方程的碰撞模型。

图2. ATA的适用范围，与量子Gibbs采样器（QGS）[12,17–20]和基于系统非久塞林德方程的非马尔可夫碰撞模型（SL算法）[39–41,46]的关系；详见第I和VA节。

我们预期上述ATA的优势使其成为一种有前景的方法，适用于涉及平衡和非平衡量子系统及非马尔可夫动力学的广泛问题的热态制备和量子模拟。此外，我们预期其对抗有限错误率的潜在鲁棒性使其成为在近中期量子硬件上实现的理想候选。

本文其余部分结构如下：在第II节中，我们介绍所研究的问题并提出第一个主要结果：高斯环境的辅助量子比特列车表示。在第III节中，我们利用该结果提出第二个主要结果：基于ATR的辅助量子比特列车量子算法。第V节为讨论部分，包括对ATA、QGS和SL算法之间关系的更详细分析。详细的误差界、估计及其证明见附录。

II. 高斯环境的辅助量子比特列车表示

在此，我们介绍可以用ATA模拟的开放量子系统模型类别（第II A节），并给出第一个主要结果：任何高斯环境都可以表示为一列辅助量子比特，依次与系统作用（第II B节）。然后我们证明该列车可以截断为有限数量的量子比特，代价是一个可控的误差，该误差在量子比特数趋于无穷的极限下消失。这一结果是ATA的基础。作为特例，我们最后证明具有高斯关联的经典噪声信号可以表示为由浴功率谱密度（PSD）的跃迁关联子（定义为浴PSD平方根的反傅里叶变换）卷积的离散二进制噪声信号。

A. 问题引入

我们考虑一个任意量子系统S，与一个或多个高斯环境耦合，以下称为浴B。高斯环境的定义见下文和附录A。我们允许系统哈密顿量$H_S(t)$是时依赖的。我们所考虑的联合系统哈密顿量为

(1)

其中$H_S$是S的哈密顿量，系统-浴相互作用$H_{SB}$包含$H$中所有同时对S和B有非平凡作用的项。一般情况下，$H_{SB}$可写为

(2)

其中$A_k$是系统的厄米算符（观测量），$B_k$是浴的厄米算符，而$g$是参数化系统-浴耦合的任意能量尺度。我们将求和中的每一项称为一个噪声通道。我们假设$B_k$是有界的。不失一般性，我们可以设$g$[50]并对每个$B_k$假设$\langle B_k \rangle = 0$；任何非零期望值可以通过修改$H_S$来消除。为简化起见，我们的讨论主要考虑上述求和仅含一项的情况，因此以下省略下标$k$。在第

**III C 一般情形的结果。** 在本节中，我们呈现一般情形的结果。此外，我们假设系统初始时刻与环境解耦，即在某一初始时刻$t=0$时系统-环境态可分离，

$$\rho_{SB}(0)=\rho_S(0)\otimes\rho_B(0), \tag{3}$$

我们的目标是获得系统的完整演化，即由其约化密度矩阵

$$\rho_S(t)=\mathrm{Tr}_B\left[\rho_{SB}(t)\right], \tag{4}$$

完整描述，其中$\rho_{SB}(t)$表示由$\rho_{SB}(0)$经哈密顿量$H$产生的时间演化，$\mathrm{Tr}_B$表示对环境Hilbert空间的部分迹。$\rho_B$的高斯性质意味着$\rho_B$关于$\hat{\Omega}$和$\hat{\Omega}^\dagger$的关联函数完全由二点关联函数确定，并按照Wick定理[44]进行结构化；详见附录A。这一性质在$\hat{\Omega}$由玻色模式的集合或连续统构成时出现，其中$\hat{\Omega}$为模式产生和湮灭算符的线性组合，且$\rho_B$关于这些模式为高斯态[44]。更广泛地说，当环境由许多有效解耦的子系统组成时，无论其内部结构如何，高斯性都会出现。因此，高斯环境无处不在，广泛出现于电磁、声子、电子、振子电子以及自旋环境等多种设置中[44,51]。高斯环境的普适性意味着它们可以用许多等价方式来表示，例如用离散耗散玻色模式的构型来表示[31,36]。在下一节中，我们将确定一种特别适合量子模拟的高效表示高斯环境的方式。

**B. 辅助粒子列（Ancilla-Train）表示。** 我们现在展示环境如何表示为一列辅助量子比特，这些量子比特依次与系统相互作用。我们首先证明环境可表示为一列与系统依次相互作用的谐振子，所有谐振子均初始化为基态。一种玻色构造量子碰撞模型的方法，类似于我们的方法，最近已在文献[43]中提出，针对具有硬紫外截止的高斯环境。随后，我们证明激发到第二激发态的概率可以忽略不计，这意味着谐振子可以被初始化为$|0\rangle$态的量子比特所替代。

首先，我们相对于$\hat{H}_S$转换到相互作用绘景[52]，引入$\hat{\Omega}_I(t)=e^{i\hat{H}_S t}\hat{\Omega}e^{-i\hat{H}_S t}$。注意$\hat{\Omega}_S$的约化密度矩阵在相互作用绘景和Schr?dinger绘景中是相同的。在相互作用绘景中，当限制在式(2)中的单项时，组合系统的哈密顿量变为

$$\hat{H}_{SB}(t)=-\hat{\Omega}_I(t)\hat{\Omega}^\dagger-\hat{\Omega}\hat{\Omega}_I^\dagger(t). \tag{5}$$

可观测量$\hat{\Omega}_I(t)$可视为经典噪声信号的量子推广[51]；我们此后将其称为环境产生的量子噪声信号。

接下来，我们利用$\hat{\Omega}_I(t)$的高斯性。这意味着在状态$\rho_B$下，$\hat{\Omega}_I$的所有关联函数按照Wick定理[44]由二点关联函数确定。此处二点关联函数由下式给出：

$$G(\tau)=\langle\hat{\Omega}_I^\dagger(t)\hat{\Omega}_I(t+\tau)\rangle_{\rho_B}. \tag{6}$$

上述事实意味着系统演化$\rho_S(t)$可以通过将式(5)中的量子噪声信号$\hat{\Omega}_I(t)$替换为任何关联函数按照Wick定理结构化且产生二点关联函数$G$的可观测量$\hat{\Omega}'_I(t)$来复现。这定义了一个广泛的可观测量等价类，它们模拟与$\hat{\Omega}_I$相同的环境。

我们现在识别两个可观测量，它们在上述意义上等价于$\hat{\Omega}_I$，并且为量子模拟提供了环境的高效表示。第一个，$\hat{\Omega}_{\text{osc}}$，是玻色模式阵列（或列）中产生和湮灭算符的含时线性组合。我们利用它来获得第二个，$\hat{\Omega}_{\text{qubit}}$，后者类似地由量子比特列构造而成，并构成了下一节中辅助粒子列方法（ATA）的基础。

我们利用文献[39]中的跳跃关联子$g(\tau)$来构造$\hat{\Omega}_{\text{osc}}$和$\hat{\Omega}_{\text{qubit}}$，其定义为

$$g(\tau)=\int_0^\infty\frac{d\omega}{2\pi}\sqrt{S(\omega)}e^{-i\omega\tau}, \tag{7}$$

其中$S(\omega)$是$\hat{\Omega}$的功率谱密度。注意$S(\omega)$始终非负，因此平方根也是非负的。跳跃关联子$g$是$G$的唯一半正定卷积平方根，即$g*g(\tau)=G(\tau)$，其中$*$表示卷积[53]：

$$G(\tau)=(g*g)(\tau). \tag{8}$$

我们通过引入一个辅助系统来构造$\hat{\Omega}_{\text{osc}}$，该系统由无限多个独立的谐振子（或自由玻色模式）$\hat{a}_k$组成，其玻色湮灭算符为$\hat{a}_k$。我们将辅助系统初始化为所有辅助模式的真空态$\rho_{\text{osc}}(0)=|0\rangle\langle 0|^{\otimes\infty}$，其中$|0\rangle$为模式$\hat{a}_k$的真空态。我们将$\hat{\Omega}_{\text{osc}}$构造为产生和湮灭算符的滑动线性组合：

$$\hat{\Omega}_{\text{osc}}(t)=\sum_{k\in\mathbb{N}}\left[g(k\Delta)^\dagger\hat{a}_{k}^\dagger+g(k\Delta)\hat{a}_{k}\right], \tag{9}$$

其中$g^\dagger(\tau)=g(-\tau)$，$\Delta$为我们选定的时间尺度。可以验证，在状态$\rho_{\text{osc}}$下，$\hat{\Omega}_{\text{osc}}$满足Wick定理，且具有二点关联函数

$$\langle\hat{\Omega}_{\text{osc}}^\dagger(t)\hat{\Omega}_{\text{osc}}(t+\tau)\rangle_{\rho_{\text{osc}}}=\sum_{k\in\mathbb{N}}|g(k\Delta)|^2+2\text{Re}\left[\sum_{m\neq n}g(m\Delta)^*g(n\Delta)e^{-i\omega(t+\tau)}\right]\Bigg|_{\omega\to(\tau)/\Delta}.$$

更准确地说，

$$\langle\hat{\Omega}_{\text{osc}}^\dagger(t)\hat{\Omega}_{\text{osc}}(t+\tau)\rangle_{\rho_{\text{osc}}}=\sum_{k\in\mathbb{N}}g(k\Delta)^*g((k+m)\Delta),$$

其中$m=\tau/\Delta$。我们将其与式(8)右侧（取离散积分）进行识别。因此，上述期望值在$\Delta\to0$极限下收敛于$g(\tau)$。由于$\hat{\Omega}_{\text{osc}}$的关联函数满足Wick定理且具有二点关联函数$G$，$\hat{\Omega}_{\text{osc}}$在$\Delta\to0$极限下复现与$\hat{\Omega}_I$相同的演化。在附录B 2中，我们对此事实提供了严格证明，并给出了选择非零$\Delta$所产生的误差界。

我们现在识别第二个可观测量$\hat{\Omega}_{\text{qubit}}$，它将环境用量子比特表示。首先注意，在$\Delta\to0$极限下，$\hat{\Omega}_{\text{osc}}$中任何辅助谐振子泄漏到第二激发态的情况都可忽略。即使$\Delta$不为零，泄漏误差也可以控制（见附录B 3）。这使我们能够在极小的误差代价下将辅助系统投影到每个谐振子的基态和第一激发态子空间中。这有效地将$\hat{\Omega}_{\text{osc}}$中的谐振子替换为量子比特。基于这一结果，我们引入量子比特辅助系统$\hat{Q}$，由无限多个量子比特组成，其Pauli算符为$\hat{\sigma}_k$。我们将辅助系统初始化为所有量子比特的$|0\rangle$态，并引入能量截断的辅助粒子列算符

$$\hat{\Omega}_{\text{qubit}}(t)=\sum_{k\in\mathbb{N}}\left[g(k\Delta)^*\hat{\sigma}_k^-+g(k\Delta)\hat{\sigma}_k^+\right], \tag{11}$$

其中$\hat{\sigma}_k^\pm=(\hat{\sigma}_x^{(k)}\mp i\hat{\sigma}_y^{(k)})/2$是量子比特$\hat{Q}_k$的常用阶梯算符。根据上述论证，在$\Delta\to0$极限下，该信号满足Wick定理且具有正确的二点关联函数$G(\tau)$。因此，$\hat{\Omega}_{\text{qubit}}$和$\hat{\Omega}_{\text{osc}}$在此极限下产生相同的系统动力学。

上述结果表明，任何高斯环境（或量子噪声信号$\hat{\Omega}_I$）都可以表示为一列量子比特，它们以跳跃关联子控制的含时方式依次与系统相互作用。量子噪声信号的辅助粒子列表示（ATR）在$\Delta\to0$极限下收敛于精确的量子噪声信号。

**经典噪声的辅助粒子列表示。** 当环境功率谱密度为偶函数时，辅助粒子列表示会出现一个具有启发性意义的特例，或者说等价地，关联函数为实值[51]。在这种情况下，$\hat{\Omega}_I$可以表示为作用于系统上的经典标量噪声信号，具有高斯关联——或称有色噪声信号——通过$\hat{\Omega}$与系统耦合。在这种情况下，$\hat{\Omega}_{\text{osc}}$为实值，式(15)变为

$$\hat{\Omega}_{\text{qubit}}(t)=\sum_{k\in\mathbb{N}}g(k\Delta)\hat{\sigma}_x^{(k)}.$$

注意$\hat{\sigma}_z^{(k)}$现在是总时间演化的运动积分。由于每个辅助量子比特都初始化为$\hat{\sigma}_z^{(k)}$两个本征态的等权叠加，对辅助系统进行部分迹等价于以各$\frac{1}{2}$的概率将$\hat{\sigma}_x^{(k)}$替换为1和$-1$。因此，我们可以将每个$\hat{\sigma}_x^{(k)}$替换为从$\{-1,+1\}$中均匀采样的独立随机变量$\epsilon_k$。这将$\hat{\Omega}_{\text{qubit}}(t)$替换为以下随机（标量）函数：

$$\Omega_{\text{classical}}(t)=\sum_{k\in\mathbb{N}}g(k\Delta)\epsilon_k. \tag{12}$$

由于上述右侧仅是一个标量随机函数，我们将它识别为经典噪声信号。辅助粒子列表示的这一特殊情形表明，经典高斯噪声信号总是可以通过将其跳跃关联子与白噪声信号进行卷积来生成。我们在图3中提供了该结果的说明。

**图3. 经典有色噪声的辅助粒子列表示：** 具有偶数功率谱密度的环境等价于作用于系统的经典噪声信号[44]。在此我们的结果表明，任何具有高斯关联的经典噪声信号——或有色噪声信号——都可以通过对一系列相同脉冲求和来生成，这些脉冲具有随机符号（红色和白色圆点）和间隔$\Delta$，并取极限$\Delta\to0$。此处脉冲由$g(t)$给出[45]。

上述结果表明，辅助粒子列表示可视为信号处理中已知结果的量子噪声推广，即有色噪声信号可以通过滤波白噪声信号来生成[45]；在此我们将跳跃关联子识别为滤波核。

**C. 将ATR截断为有限辅助粒子数。** 我们现在证明，ATR中的无限辅助粒子列可以被具有有限寄存器大小的辅助粒子列量子噪声信号很好地近似：我们证明量子噪声信号$\hat{\Omega}_I$可以被模拟噪声信号替代，其形式如式(15)所示，通过将求和截断为有限项，即选取有限步长$\Delta$并丢弃求和中$k\Delta>T$的项（其中$T$为某截止时刻）来实现。

我们在附录B 2中证明，ATR可以按照上述方式截断为有限辅助寄存器，代价是任意小的误差，前提是我们同时将$g$替换为滤波后的跳跃关联子[54]，定义为

$$g_{\text{filt}}(\tau)=\int_{-\infty}^{\infty}f(\omega)g(\tau)e^{i\omega\tau}d\omega, \tag{13}$$

其中$f(\omega)$为低通滤波核，满足$f(\omega)=1$（对$|\omega|\leq\omega_c$）和$f(\omega)=0$（对$|\omega|\geq\omega_{c'}$）。在这两个区间之间，$f$应几乎处处具有均匀有界的导数和二阶导数，这意味着$g_{\text{filt}}$的衰减尺度大约与$g$相同——在附录B 2中我们通过多项式插值提供了$f$的一个具体构造。以下，我们给出三个算符，当选择足够小的$\Delta$和足够大的截止时刻$T$时，它们能够忠实地复现$g_{\text{filt}}$产生的演化：

$$\hat{\Omega}_{\text{osc,filt}}(t)=\sum_{k\in\mathbb{N}}\left[g_{\text{filt}}(k\Delta)^*\hat{a}_k^\dagger+g_{\text{filt}}(k\Delta)\hat{a}_k\right], \tag{14}$$

$$\hat{\Omega}_{\text{qubit,filt}}(t)=\sum_{k=0}^{N_{\text{max}}}\left[g_{\text{filt}}(k\Delta)^*\hat{\sigma}_k^-+g_{\text{filt}}(k\Delta)\hat{\sigma}_k^+\right], \tag{15}$$

$$\hat{\Omega}_{\text{ATA}}(t)=\sum_{k=0}^{N_{\text{max}}}\left[g_{\text{filt}}(k\Delta)^*\hat{\sigma}_k^-+g_{\text{filt}}(k\Delta)\hat{\sigma}_k^+\right]\mathbb{1}_{|k-t/\Delta|\leq N_{\text{win}}}, \tag{16}$$

其中$N_{\text{max}}=\lfloor T/\Delta\rfloor$。我们将使用前两个可观测量作为中间步骤来界定用$\hat{\Omega}_{\text{ATA}}$替代$\hat{\Omega}_I$的误差，后者只需要访问有限大小为$N_{\text{win}}$的辅助寄存器，并作为辅助粒子列方法（ATA）的基础。以下，我们将此误差定义为由具有正$\Delta$和有限$T$的量子噪声信号$\hat{\Omega}_{\text{ATA}}$生成的演化$\rho_S^{\text{ATA}}$与由$\hat{\Omega}_I$产生的真实演化$\rho_S^{\text{true}}$之间的迹范数距离。

我们的误差界依赖于环境的两个量：特征系统-环境相互作用速率$\gamma$和环境特征关联时间$\tau_B$，定义为

$$\gamma=\left\|\hat{\Omega}_I(t)\right\|, \quad \tau_B=\int_0^\infty|G(\tau)|d\tau, \tag{17}$$

其中$\|\cdot\|$表示算符范数。重要的是，$\gamma$定义了环境诱导演化速率的上界[39]。

以下，我们确定$\Delta$和$T$的值，使得截断引入的误差小于某个给定的目标相对精度$\epsilon$。具体来说，我们要求

$$\|\rho_S^{\text{ATA}}-\rho_S^{\text{true}}\|_1\leq\epsilon\|\rho_S^{\text{true}}\|_1, \tag{18}$$

其中$\|\cdot\|_1$表示迹范数。我们采用三角不等式将此误差界定为三个贡献之和：第一个误差$\epsilon_{\text{filt}}$来自用$g_{\text{filt}}$替代$g$。第二个误差$\epsilon_{\text{qubit}}$来自用$\hat{\Omega}_{\text{qubit,filt}}$替代振荡子信号$\hat{\Omega}_{\text{osc,filt}}$。最后一个误差$\epsilon_{\text{tail}}$来自丢弃涉及超过截止时刻$T$的项，得到$\hat{\Omega}_{\text{ATA}}$。我们在以下分别讨论这三个误差。

离散化误差$\epsilon_{\text{filt}}$可以由依赖于环境功率谱密度有效紫外截止的量来界定，即在某个能量以上$S(\omega)$有效消失。具体来说，我们在附录B 2中发现，选择$\Delta$使得

$$\Delta\leq\frac{\epsilon^2}{C\gamma\omega_c}, \tag{19}$$

时用$g_{\text{filt}}$替代$g$对$\rho_S$引入的修正的迹范数阶数最多为$C\gamma\omega_c\Delta$。特别是，当$\omega_c$和$\gamma$可积时，$\epsilon_{\text{filt}}\propto\omega_c\Delta$。对于这样的功率谱密度，我们通过选择以下紫外截止来达到目标精度：

$$\omega_c\geq\frac{\epsilon}{\gamma\tau_B}, \tag{20}$$

确保$\omega_c$的要求是为了保证频率滤波不会显著改变$\hat{\Omega}_I$和$G$（详见附录C）。我们期望$\omega_c$可以由与耦合到环境的可观测量有效谱宽来选择，这等价于算法中的Trotter步长，因此$\omega_c$可以选择与Trotter步长同阶。

我们现在考虑用量子比特替代谐振子产生的误差$\epsilon_{\text{qubit}}$。当谐振子泄漏到第二激发态可以忽略时，此误差很小。当$\Delta\to0$时泄漏消失，因为与每个辅助量子比特相互作用的强度变得可以忽略。在附录B 3中，我们将这一直觉精确化并证明当选择$\Delta\leq\epsilon/(\gamma T)$时误差保持在目标精度范围内。

最后，我们考虑丢弃距离给定时刻$t$超过截止$T$的辅助量子比特产生的误差$\epsilon_{\text{tail}}$。在附录B 4中，我们证明当选择$T\geq\frac{1}{\gamma}\log\frac{1}{\epsilon}$时$\epsilon_{\text{tail}}$在目标精度范围内。

综合以上结果，当选择

$$\Delta\leq\frac{\epsilon^2}{C\gamma\omega_c}, \quad \omega_c\geq\frac{\epsilon}{\gamma\tau_B}, \quad T\geq\frac{1}{\gamma}\log\frac{1}{\epsilon} \tag{21}$$

时，$\hat{\Omega}_{\text{ATA}}$以相对目标精度$\epsilon$复现正确的动力学。在这些参数下，在任一给定时刻$\hat{\Omega}_{\text{ATA}}$仅涉及有限数量的辅助量子比特，其数量为

$$N_{\text{win}}\sim\frac{T}{\Delta}\sim\frac{C\gamma\omega_c}{\epsilon^2}\cdot\frac{1}{\gamma}\log\frac{1}{\epsilon}. \tag{22}$$

重要的是，我们注意到当环境功率谱密度衰减足够快时，所需辅助量子比特数随$\epsilon$呈低次多项式增长。

上述结果表明，任何高斯环境都可以被$\hat{\Omega}_{\text{ATA}}$准确地表示，即使在具有无限多个无界自由度的情况下也是如此。这具有重要意义，因为$\hat{\Omega}_{\text{ATA}}$是由有限辅助寄存器中的$N_{\text{win}}$个量子比特构造而成的。

重要的是，有界；对于上述的选择，[55]。这一事实使得利用辅助粒子-列车表示对系统演化进行高效的 Trotter 化成为可能。在下一节中，我们将利用这一事实构造一个高效的量子算法，用于模拟与任意高斯环境在任意耦合强度下相连接的开放量子系统的动力学。

III. 辅助粒子-列车算法

我们现在利用 ATR 来获得我们的主要结果：一个量子算法，用于系统地模拟与任意配置的高斯环境在任意耦合强度下相耦合的开放量子系统。ATA 是通过对 ATR 中联合系统-辅助粒子（）演化进行 Trotter 化而获得的，其中使用了截断的辅助粒子-列车量子噪声信号。我们在图 4 中提供了 ATA 的紧凑伪代码。

**图 4. 辅助粒子-列车算法的伪代码。** 输入参数包括初始系统状态、热浴功率谱密度、系统哈密顿量、系统耦合算符、系统-热浴耦合、总模拟时间以及（相对）目标精度。由（26）给出的表示由系统哈密顿量生成的 Trotter 相互作用，由（27）和（29）给出的表示由辅助粒子-系统耦合生成的时间演化中的 Trotter 迭代。我们用表示的期望值为。

本节结构如下：为简便起见，我们首先介绍单噪声通道的 ATA（III A 节）[56] 并提供资源估计（III B 节）。然后在 III C 节中将结果推广到高斯环境的最一般（多通道）配置。

**A. ATA 的构造**

ATA 模拟 II A 节中定义的一类开放量子系统的演化：即任意物理系统，具有（可能含时哈密顿量），通过任意系统-热浴耦合连接至任意配置的高斯环境。为简便起见，我们首先考虑可分解为系统观测量与热浴观测量乘积的情况：，其中。这里参数化系统-热浴耦合，可以任意大。我们在 III C 节中考虑更一般的情况。

ATA 在系统量子比特寄存器上操作，编码系统状态，以及辅助粒子寄存器，用于通过 ATR 模拟热浴，如 II 节所述。系统可以通过标准的哈密顿量模拟技术编码在 [2,57–62]。此后，我们简单地用替代，而和指代表示其在被模拟系统中的对应物。

ATA 以在模拟中的初始状态 [63] 作为输入。然后 ATA 通过迭代地通过 Trotter 化演化将系统沿时间向前传播来模拟由引起的演化。具体来说，ATA 的第次迭代返回的一个状态，其对应的约化密度矩阵

(23)

提供了的控制近似（见下文）。这里和，称为 Trotter 步，是模拟的输入参数，控制精度和输出分辨率。

ATA 的第一步将完整的系统-辅助粒子寄存器准备在初始状态，其中是中的个状态的张量积。然后 ATA 通过与迭代地生成轨迹。具体来说，通过一系列门序列——Trotter 酉算符——迭代获得，其中第个 Trotter 酉算符，将系统从时间向前传播到：

(24)

我们使用二阶 Trotter 化方案，使得由

(25)

给出，其中和是由和生成的 Trotter 酉算符：

(26)

(27)

其中表示时间序符号，表示来自式（15）的模拟热浴信号

(28)

这里，其中和是模拟的输入参数，与一起控制模拟精度；更多细节见下文。

Trotter 酉算符的第一个门序列可以通过现有的哈密顿量模拟技术实现 [2,57–62]。热浴传播子本身是由对易酉算符的乘积构成的，每个都是由将量子比特中的单个量子比特与耦合生成的：其中

(29)

其中；表示积分区域的边界，同时考虑了有限所带来的限制以及与 Trotter 步对应的时间区间的界限。注意，对于局域的或准局域的，可以有效地编译为有限量子门序列，甚至可能是经典可处理的。因此，我们期望的编译通常成本低廉。

在附录 B 中，我们表明当和选择得足够小时，ATA 对任意非零目标精度返回真实的系统演化。我们通过结合 II 节的讨论与时间演化二阶 Trotter 化的误差界 [64] 来建立这一结果。我们的界假设系统哈密顿量不含时，[65]，但我们期望推广到含时哈密顿量是直接的。具体来说，我们表明 ATA 的输出轨迹保证能以给定的（相对）目标精度重现真实的系统演化，

(30)

条件为

(31)

其中，而，和定义于式（17）和（20）。对于非平凡情况的更详细的约束条件见附录 B 5。结果通过结合 II 节中的误差界与演化的 Trotter 化的误差界得到；细节见附录 B。

我们最后强调，中的辅助量子比特只需在它们被使用的第一个 Trotter 步中准备：即辅助量子比特只需在 Trotter 步准备。同样，辅助量子比特在其最后一次影响系统演化的 Trotter 步之后可以丢弃。取决于硬件，特别是量子比特的资源成本，丢弃的辅助量子比特可以被重置并重新引入寄存器。在这种方案中，只需持有量子比特；在这些中，每个 Trotter 步中有被重置并重新引入寄存器，而其余量子比特被保留。仍然可以在硬件中实现辅助粒子寄存器作为持续传入量子比特的流，在这种情况下资源效率更高 [66,67]。

**B. 资源成本**

我们现在估计 ATA 与由生成的标准酉时间演化相比的额外资源成本。

我们首先考虑实现 ATA 所需的额外量子比特数——即辅助粒子寄存器中的最少量子比特数。回顾式（31）下的讨论，ATA 只需访问辅助粒子寄存器中的量子比特。就误差缩放而言，所需辅助量子比特数上界为

在极限下，所需辅助量子比特数的缩放为

(32)

接下来，我们考虑实现 ATA 与标准酉演化相比的门复杂度。每个 Trotter 步涉及个酉算符的乘积。如 III A 节所述，每个由与量子比特中可观测量耦合生成。因此，我们期望可以有效地编译为有限量子门序列。与标准 Trotter 化哈密顿量模拟相比，ATA 的额外门计数缩放为，由此导致

(33)

如果我们选择与成比例，即对整数，那么门与相同，只是辅助粒子寄存器有一个位移，其中" "表示最小公倍数。在这种情况下，门可以被高效地预编译。给定一组用于模拟孤立系统演化的酉门，的编译很容易跟随。对于不含时的，所有相同，因此预编译的数量与模拟时间无关。

我们在式（32）和（33）中对 ATA 的资源估计表明，当热浴功率谱密度衰减足够快时，相对于标准 Trotter 化演化的额外资源成本是目标精度倒数的低次多项式。重要的是，额外成本与哈密顿量模拟问题的复杂度无关。因此，对于广泛的问题类别，用 ATA 替代标准 Trotter 化演化的额外成本可以忽略不计。

**C. 多通道环境的定义**

我们最后展示如何将 ATA 直接扩展到任意高斯环境和任意系统-热浴耦合。

不失一般性 [39]，在这种情况下，我们可以将表示为

(34)

这里和分别代表列向量和，其中表示和上的观测量，表示噪声通道的总数。我们使用这种向量表示法来强调与单通道情况下的相似性。热浴关联函数现在变为矩阵值的，，其中元素由

(35)

给出

我们将跳跃关联子的定义扩展为的正类型的唯一卷积矩阵平方根，即，其中表示沿负实轴有割线的矩阵平方根 [68]。

有了这些定义，我们将 ATA 从单通道扩展到多通道情况如下：首先，将辅助粒子寄存器从一条扩展为条辅助粒子列车，即具有泡利算符的量子比特数组 [69]。从这些辅助粒子列车中，我们可以用截断的模拟噪声信号替代每个，其中

(36)

其中表示滤波后跳跃关联子的第个条目。注意，每条辅助粒子列车不一定模拟单个通道；相反，每个通道通常由所有辅助粒子列车的组合生成。

多通道情况的 ATA 定义与单通道情况相同，只需将系统-辅助粒子耦合替换为考虑多个通道的形式：

(37)

该酉算符可以通过利用不同辅助粒子之间的对易性进一步 Trotter 化（见附录 D）；如果相互对易，则直接分解为作用于单个辅助粒子的算符。

除了上述替换外，多通道情况的 ATA 与单通道情况不变。事实上，我们在附录 B 中的误差界是为多通道情况推导的；它表明如果，且满足式（31），则多通道热浴的 ATA 保证以式（30）的意义下具有精度，其中，和现在分别通过将，和用，和替换来定义。这里，对于任何矩阵，表示的矩阵范数。

IV. 数值模拟

在本节中，我们模拟了 ATA 用于一个自旋-玻色子模型，该模型由一个双能级系统与具有欧姆谱密度的环境耦合组成。我们的目标是展示算法的工作原理并探究其资源需求，在实际中这些需求可能比 III B 节中的严格界小得多，后者覆盖了最坏情况。

我们模拟的模型哈密顿量为

(38)

其中是具有欧姆截止的高斯环境观测量。功率谱密度为

(39)

其中是热浴的逆温度，是高斯截止。在所有模拟中，我们设置，和。我们考虑三种不同温度：，和=10。

我们考虑通过将系统初始化为所产生的动力学。通过直接时间演化，我们使用多达 13 个辅助量子比特模拟 ATA 下的系统演化；伪代码见图 5。在所有模拟中，我们设置以最小化数值成本 [70]。从 III B 节中和式（31）中的表达式，我们得到以下不等式：

(40)

其中是截止时间和是紫外截止。在所有模拟中，我们选择的最大可能值。图 5 的前两列展示了和的选择方式，使得跳跃关联子和功率谱密度在各自截止范围外的权重很小。在所有模拟中，我们根据辅助量子比特数和在截止时间的约束下，选择最小的可能辅助量子比特间距，，通过。这里是一个小的无理数修正，用于避免 Trotter 化演化中的 Floquet 共振 [71]。为了数值可处理性，我们最后将辅助粒子寄存器的希尔伯特空间限制为最多四个同时激发的子空间。

**图 5. 欧姆自旋-玻色子模型的 ATA 模拟。**（细节见正文）。在所有图中，我们设置，和。三行分别对应不同的热浴逆温度：（上）、（中）和（下）。第 1 列：跳跃关联子的绝对值。灰色条显示模拟中所选截止时间。第 2 列：所选参数的热浴功率谱密度。灰色条显示不同辅助量子比特数下的紫外截止。最内层区域条对应 5 个辅助量子比特，最外层区域对应 13 个辅助量子比特。

第三列：使用不同数量辅助量子比特进行模拟时，自旋处于态的布居数随时间的演化；标记使用个辅助量子比特的模拟。第四列：使用13个辅助量子比特的模拟与使用较少辅助量子比特的模拟之间的迹距离。我们看到，当使用13个辅助量子比特时，各模拟相互收敛至2%以内。在图5的第三列中，我们绘制了ATA模拟中态布居数的演化。在第四列中，我们绘制了系统态演化与使用13个辅助量子比特的演化之间的迹距离。这作为关于辅助量子比特数收敛性的指标。如图5所示，当比较使用11和13个辅助量子比特的模拟时，ATA模拟收敛至以内。可以清楚地看到，随着辅助量子比特数增加，距离减小——这一趋势在瞬态过程中最为明显。这表明算法关于辅助量子比特数是收敛的，对于13个辅助量子比特的模拟，有限辅助量子比特数带来的误差为百分之几。将我们的结果与第III B节的误差界进行比较是有启发性的。由式(32)和式(30)，我们看到确保总误差小于所需的辅助量子比特数为(41)，其中是总模拟时间。对于我们探测的参数，我们观察到模型在时间内弛豫到稳态；见图5。为了确保在此时间跨度内误差低于，我们的解析界表明我们需要个辅助量子比特，如表II第三列所示。相比之下，我们的模拟表明，对于个辅助量子比特，收敛在百分之几以内。因此，这些结果表明，实际上，数值收敛所需的辅助量子比特数远少于我们的界所预测的。表II. 欧姆自旋-玻色子模型的参数和辅助量子比特数的解析上界。该估计基于在处的10%误差阈值。辅助量子比特界：0.1，0.2，2.7，0.5%；1，0.2，1.2，1.6%；10，0.7，1.9，1.2%。V. 讨论。在本文中，我们介绍了一种量子算法，用于系统地模拟与任意高斯环境以任何耦合强度相互作用的量子系统：辅助量子比特列车算法（ATA）。ATA可用于模拟平衡和非平衡问题的动力学和稳态，并且适用于非马尔可夫区域。我们提供了算法精度的严格界，并表明通过适当选择模拟分辨率参数，模拟误差可以被任意地减小。ATA的资源成本是关于逆目标精度的多项式，前提是浴功率谱密度（PSD）具有适度的衰减特性，并且由模拟模型的四个物理能量尺度的组合控制：与环境耦合的系统可观测量在海森堡图像中的特征速度，以及浴关联时间、系统-浴耦合强度和环境的有效紫外截断。重要的是，资源成本与系统的复杂性无关。因此，对于广泛的问题，与模拟不存在外部环境时的时间演化的成本相比，ATA只需要少量的额外资源。量子噪声的辅助量子比特列车表示（ATR）提供了我们工作的一个底层结果：任何高斯量子环境都可以通过将量子二进制白噪声信号——即辅助量子比特列车——与一个时间局域卷积核进行卷积来模拟。一个推论是，具有高斯关联的经典噪声信号可以通过将白噪声信号与跳跃关联函数进行卷积来生成。这一结果与经典信号处理中的已知结果一致，即有色噪声信号可以通过对白噪声信号进行滤波获得[45]；这里滤波核由跳跃关联函数给出。因此，ATR通过非马尔可夫碰撞模型将此结果推广到量子噪声。对于量子噪声，滤波核（跳跃关联函数）简单地从实值函数推广为复值函数，而经典（二进制）白噪声信号——比特流——则被替换为量子比特流。A. 与量子吉布斯采样和系统性Lindblad方程算法的关系。我们现在讨论ATA与文献中最近出现的一些开放量子系统量子模拟算法之间的关系：QGS[12,16–20]和基于系统性严格推导的Lindblad方程的碰撞模型（SL）[39–41]。QGS算法给出给定模型哈密顿量的精确吉布斯态。这是通过基于具有准局域跳跃算符且表现出KMS细致平衡的Lindblad过程的碰撞模型实现的[46][72]。已知与热平衡环境相连的物理系统在弱系统-浴耦合极限下会收敛到这样的精确吉布斯态。然而，对于非零有限耦合——这在物理系统中是存在的——由于为分区量子系统定义适当的子系统哈密顿量的模糊性，稳态与吉布斯态之间存在由系统-浴耦合控制的修正偏差[48]。因此，从量子模拟的角度来看，QGS算法仅在系统-浴耦合消失的极限下才能给出开放量子系统的准确物理稳态。与QGS算法的出现并行的是，一系列系统性严格推导的Lindblad方程的出现，它们绕过了旋转波（或久期）近似。这些Lindblad方程已被成功应用于耗散量子多体系统的经典模拟中[73–78]。与QGS Lindblad过程类似，这些系统性Lindblad方程具有准局域跳跃算符，因此可以使用碰撞模型类型的方法高效地实现为量子算法[46]。我们将此类量子算法称为SL算法。与QGS算法不同，SL算法从底层微观模型系统地推导而来，并扩展到非平衡问题，例如涉及驱动或非平衡储库的问题。这些SL方法给出的动力学和稳态——平衡或非平衡——精确到与系统-环境耦合成正比的修正（相对于环境的内在关联时间尺度[49]）。因此，从量子模拟的角度来看，SL算法给出的稳态在相同的近似水平上准确——即QGS和SL算法都在系统-浴耦合消失的极限下是准确的。ATA补充了上述QGS和SL算法，提供了重要的优势：(1) ATA模拟的精度仅受可用资源的限制；它可以以任意精度给出模拟微观模型的精确稳态和动力学。(2) ATA仅涉及简单的严格局域系统-辅助量子比特耦合，因此不一定需要编译准局域系统-辅助量子比特耦合。特别是，如果系统中的耦合算符可以表示为严格局域算符，则酉系统-辅助量子比特耦合算符也是严格局域的。这与许多基于马尔可夫碰撞模型的量子算法形成对比[17,39]。(3) ATA不需要编译准局域Lamb位移或其推广。对于给定的固定目标精度，ATA的额外成本是一个大小与系统-浴耦合线性标度的辅助量子比特寄存器。因此，当电路深度是限制资源时，ATA在适用范围更广和实现更简单方面都优于QGS和SL算法。在本文准备过程中，类辅助量子比特列车的构造出现在QGS文献中[19–21]。这里的目标是模拟一个以热（吉布斯）态作为近似稳态的Lindblad过程。参考文献[19–21]为该构造提供了严格的性能保证，证明了辅助量子比特列车广泛应用于QGS，并且在其中带来了严格局域系统-辅助量子比特耦合的优势。B. 噪声鲁棒性。ATA的一个最终潜在优势，也是其他开放量子系统模拟算法所共有的[12,17–20,39,40]，是其由于提取熵的本性而对硬件错误具有潜在的内在鲁棒性[79]。具体来说，执行ATA的硬件中的噪声将映射为模拟微观模型中的噪声信道。由于模拟已经涉及耗散机制，该噪声信道将与其他耗散和噪声信道竞争，因此如果足够弱，则只会产生输出（模拟稳态）的连续扰动：因此，对于诸如ATA的开放量子系统模拟器，噪声引起的误差通常是噪声强度的有限斜率函数。等价地，由于ATA模拟中的耗散信道为稳态扰动提供了有限的弛豫时间，如果错误的特征速率相对于弛豫时间足够小，则外部噪声将产生可忽略的误差。开放量子动力学模拟中固有的误差鲁棒性的程度仍然是一个开放的研究领域。然而，最近的一些理论和数值结果印证了一个直觉，即开放量子系统动力学的模拟比闭合量子系统的模拟对错误的鲁棒性显著更强[80–83]。C. 展望。总之，ATA提供了一种在任意耦合强度和任意高斯环境的非平衡或平衡配置下系统地准确模拟开放量子系统的方法。将ATA限制在高斯环境可以看作是环境多时间关联函数的二阶累积量展开。ATA向高斯环境之外的自然扩展将是包含更高阶的累积量。我们将这样的研究留给未来的工作。我们预计该算法对硬件错误具有内在的鲁棒性，且相对简单易实现。我们预计这些结果使ATA成为非平衡量子系统量子模拟的有前景的候选算法，并可在近期量子硬件上实现。