来源论文: https://arxiv.org/abs/2605.22459v1 生成时间: May 26, 2026 18:43

0. 执行摘要

在化学物理、光合作用光致激发能传递以及凝聚态量子物理领域,精确描述有限温度下与复杂环境耦合的分子体系动力学是一项极具挑战性的任务。传统的非马尔可夫动力学模拟方法(如 HEOM、QUAPI 等)通常需要针对每一种特定初始状态进行单独的、计算极其昂贵的时间传播。为了彻底解决这一计算瓶颈,Raffaele Borrelli 和 Hideaki Takahashi 在其最新工作中提出了一种全新的一体化计算框架。

该方法巧妙地融合了三大核心理论和技术支柱:蔡-雅米奥科夫斯基同构(Choi-Jamiołkowski Isomorphism)热场动力学(Thermofield Dynamics, TFD)纯化技术,以及张量列(Tensor-Train, TT)波函数传播技术。通过将系统嵌入到一个包含“系统副本(Ancilla)”的极大纠缠态中,并在双倍温度场 Hilbert 空间中进行单次幺正动力学传播,研究者能够直接重构出完整的、包含所有初始状态演化信息的有限温度约化动力学映射 $\Phi(t)$。本文将对该项工作的理论机制、技术难点、数值 Benchmark、代码实现方案以及局限性进行深度的学术解析。


1. 核心科学问题,理论基础,技术难点,方法细节

1.1 核心科学问题:非马尔可夫约化动力学的“多状态计算瓶颈”

描述开放量子系统的演化通常依赖于约化系统密度矩阵 $\rho_S(t)$ 的传播。由于实际分子(如光合捕光复合物、有机半导体)不仅与电子自由度耦合,还深度耦合到结构化的、具有长记忆时间的振动模式(即声子浴)中,其动力学具有强烈的非马尔可夫(Memory)效应。传统的非马尔可夫动力学模拟在面对不同初始条件时(如研究多波长脉冲激发、不同的相干初始态叠加),必须对每一个初始状态重复进行昂贵的计算。这在系统维度 $d_S$ 增大时面临严重的维数灾难,因为要完全重构系统的动力学映射 $\Phi(t)$,传统方法需要至少进行 $O(d_S^2)$ 次独立模拟。

如何通过单次动力学传播,获取有限温度下包含任意因式分解初始态演化信息的完整通道(Channel)表征,是该领域的长期痛点。这就是本研究力图攻克的非马尔可夫约化映射构建难题。

1.2 理论基础一:Choi-Jamiołkowski 同构

约化动力学映射(简称量子通道) $\Phi(t, t_0)$ 描述了系统从初始时刻 $t_0$ 到 $t$ 时的状态演化:

$$\rho_S(t) = \Phi(t, t_0)[\rho_S(t_0)] = \text{Tr}_B [ U(t, t_0) (\rho_S(t_0) \otimes \rho_B) U^\dagger(t, t_0) ]$$

其中 $\rho_B$ 为环境(浴)的热平衡态。根据 Choi-Jamiołkowski 同构理论,一个定义在 Hilbert 空间 $\mathcal{H}_S$ 上的量子通道 $\Phi$ 可以等价地由一个定义在联合空间 $\mathcal{H}_S \otimes \mathcal{H}_A$ 上的正算符——Choi 矩阵 $J_\Phi$ 来唯一表示。这里,辅助系统 $A$(Ancilla)是系统 $S$ 的一个等维度副本($d_A = d_S$)。

首先构造系统与辅助系统之间的极大纠缠态(无归一化):

$$|\Omega\rangle = \sum_{i=1}^{d_S} |i\rangle_S \otimes |i\rangle_A$$

对应的 Choi 矩阵 $J_\Phi(t)$ 形式为:

$$J_\Phi(t) = \sum_{i,j} \Phi(t)[|i\rangle\langle j|_S] \otimes |i\rangle\langle j|_A$$

这是一个 $d_S^2 \times d_S^2$ 的矩阵。一旦求得 $J_\Phi(t)$,对于任意给定的初始密度矩阵 $\rho_S(0)$,其在 $t$ 时刻的响应可直接通过简单的张量收缩实现,无需再次运行任何量子动力学模拟:

$$(\Phi(t)[\rho])(m, n) = \sum_{i,j} J_\Phi(m, i; n, j) \rho(i, j)$$

为了用单次动力学模拟构建 $J_\Phi(t)$,论文构造了如下扩展空间中的投影算符 $Y(t)$:

$$Y(t) = (U \otimes I_A) (|\Omega\rangle\langle \Omega| \otimes \rho_B) (U^\dagger \otimes I_A)$$

在扩展 Hilbert 空间中,对其进行环境偏迹收缩,即直接得到 Choi 矩阵:

$$J_\Phi(t) = \text{Tr}_B (Y(t))$$

1.3 理论基础二:热场动力学(TFD)有限温度纯化

在有限温度下,环境状态 $\rho_B$ 是一个混态(canonical ensemble),这使得基于波函数的传播方法无法直接使用。传统上需采用密度矩阵传播或蒙特卡洛抽样。为了保持波函数方法的高效性,论文引入了热场动力学(TFD)。

TFD 通过引入环境系统 $B$ 的一个虚拟镜像副本 $\tilde{B}$(tilde 空间),将混态 $\rho_B$ 纯化为扩展空间 $\mathcal{H}_B \otimes \mathcal{H}_{\tilde{B}}$ 中的纯态 $|\rho_B^{1/2}\rangle_{B\tilde{B}}$:

$$|\rho_B^{1/2}\rangle_{B\tilde{B}} = \sum_k \rho_{B,kk}^{1/2} |k\rangle_B \otimes |k\rangle_{\tilde{B}}$$

对于简谐声子浴($H_B = \sum_k \omega_k a_k^\dagger a_k$),纯化状态可以通过对 $T=0$ 的真空态 $|0\rangle_{B\tilde{B}}$ 施加温度相关的 Bogoliubov 变换得到:

$$|\rho_B^{1/2}\rangle_{B\tilde{B}} = e^{i G(\beta)} |0\rangle_{B\tilde{B}}$$$$G(\beta) = i \sum_k \theta_k(\beta) \left( \tilde{a}_k a_k^\dagger - a_k \tilde{a}_k^\dagger \right)$$

其中 $\tanh\theta_k(\beta) = e^{-\beta \omega_k / 2}$,$\beta = 1/(k_B T)$。

为了数值计算的便利,作者使用了一个极具智慧的佯攻:将温度依赖性从初始状态转移到哈密顿量上。通过幺正变换 $H_\theta = e^{-i G(\beta)} H e^{i G(\beta)}$,我们可以在零温真空态 $|0\rangle_{B\tilde{B}}$ 下,利用温度调制的哈密顿量 $H_\theta$ 进行实时间演化。对于线性振动电子耦合模型,变换后的温度依赖哈密顿量 $H_\theta$ 具有如下紧凑形式:

$$H_\theta = \sum_{n,m} \left( \varepsilon_{nm} + \sum_{k=1}^{N_B} \frac{g_{knm}(eta)}{\sqrt{2}} (b_k^\dagger + b_k) \right) |n\rangle\langle m| + \sum_{k=1}^{N_B} \omega_k b_k^\dagger b_k$$

这里,包含物理模式和镜像虚拟模式的算符被统一写作 $b_k$ 与 $b_k^\dagger$,其有效的温控耦合系数 $g_{knm}(eta)$ 和频率 $\omega_k$(频率可正可负,负值对应伴随的 tilde 模式)直接体现了有限温度效应。这使得整体演化态写为:

$$|\varphi(t)\rangle = e^{-i H_\theta t} |\Omega\rangle |0\rangle_{B\tilde{B}}$$

而 Choi 矩阵则由下式给出:

$$J_\Phi(t) = \text{Tr}_{B\tilde{B}} [|\varphi(t)\rangle \langle \varphi(t)|]$$

1.4 理论基础三:张量列(Tensor-Train, TT)表示与 TDVP 传播

即使经过了纯化,由于声子浴模式数 $N_B$ 极其庞大(论文中多达数百个模式),多体波函数 $|\varphi(t)\rangle$ 仍会面临维数灾难。为此,作者引入了张量列(TT)格式(在物理学中常称为矩阵乘积态, MPS)。

一个包含 $L = D + N_B + 2$ 个自由度的整体波函数被写为如下矩阵乘积形式:

$$|\varphi\rangle = \sum_{i_1, \dots, i_L} G^{(1)}(i_1) G^{(2)}(i_2) \cdots G^{(L)}(i_L) |i_1\rangle \cdots |i_L\rangle$$

其中 $G^{(k)}(i_k)$ 是三阶核心张量,其维度为 $r_{k-1} \times p_k \times r_k$,而 $r_k$ 即为 TT 秩(Rank),控制着量子态中纠缠剪裁的精度。为了最大化 TT 矩阵乘积表示的精度并控制计算开销,系统自由度的排序设计至关重要。论文采用了如下最优链式排序:

$$z_1 - \cdots - z_D - A - S - b_1 - \cdots - b_{N_B}$$
  • $z_i$ 为用于处理静态无序的辅助零频玻色模式(见下文);
  • $A$ 为辅助系统(Ancilla),$S$ 为物理主系统;
  • $b_j$ 为热声子浴物理与镜像虚拟模式,按照其频率绝对值 $|\omega_j|$ 从小到大升序排列。

时变变分原理(TDVP): 研究者利用单位置(one-site)时变变分原理(TDVP)在切空间上投影演化薛定谔方程。单位置 TDVP 在 TT 传播中具有优异的模守恒(Norm-preserving)特性,特别适合具有星形耦合拓扑(即主系统与大量独立浴模式耦合)的系统,在保证中等 TT 秩(最大 $\sim 220$)的前提下,能够实现数皮秒的高精度演化。

1.5 技术难点与创见

  1. 静态无序的一体化处理:实际体系中常伴随 site-energy 静态无序。传统方法需要通过蒙特卡洛抽样产生几百个随机哈密顿量进行重复动力学传播后平均。本工作创造性地将静态无序分布 $p(\sigma)$ 编码为具有零频率的辅助玻色模式 $z_{nm}$,将其直接写入变分哈密顿量 $\bar{H}_\theta$: $$\bar{H}_\theta = H_\theta + \sum_{n,m} V_{nm} \frac{\sigma_{nm}}{\sqrt{2}} (z_{nm}^\dagger + z_{nm})$$ 由于零频模式不参与系统和真实声子浴之间的能量交换,在初始态中将其设置为对应高斯分布的相干态,便可在单次波函数演化中自动实现系综平均,无需任何显式抽样,极大节省了计算资源。
  2. 迹守恒的非退化性:尽管变分 TDVP 传播能够严格保证波函数模守恒 $\langle\varphi(t)|\varphi(t)\rangle = 1$,但是在被投影到有限 rank 的张量流流形(TT manifold)上时,偏迹收缩后不一定完美保证约化通道的迹守恒(Trace Preservation, TP),即不严格保证 $\text{Tr}_S J_\Phi = I_A$。作者通过严格的步长测试,证实了通过精细调节演化步长 $\Delta t$(例如从 $1\text{ fs}$ 缩减至 $0.25\text{ fs}$),TP 残差能够从 $1.3 \times 10^{-4}$ 锐减至 $4.0 \times 10^{-8}$,且几乎与 TT 秩无关,这从数值上完美打消了理论学界的顾虑。

2. 关键 Benchmark 体系,计算所得数据,性能数据

2.1 Benchmark 系统:Fenna-Matthews-Olson (FMO) 复合物

作者选择光合作用领域最具标志性的 Benchmark 系统——FMO 复合物的亚基 A 进行了模拟。该系统由 8 个拟菌绿素(BChls 1-8)分子组成($d_S = 8$),是一个非常复杂的振动电子强耦合开放系统。系统哈密顿量采用 Moix 等人提出的经典参数,每个 site 都耦合了具有丰富细节的 site-specific 结构化光谱密度(Spectral Density)。

2.2 光谱密度的低秩离散化与 TFD 系统自由度

为了将连续的物理光谱密度 $J(\omega)$ 转化为一组合适的有限个离散模式,作者使用了他们之前发展的低秩插值分解(Interpolative-Decomposition, ID)方案。在温度 $T=300\text{ K}$、时间上限 $T_c = 800\text{ fs}$ 的边界条件下,每个 BChl site 的局域热噪声光谱密度 $S_\beta(\omega) = J(\omega)[1 + \coth(\beta \omega / 2)]$ 被精确地离散化为 50 个有效模式(包括物理和镜像 tilde 模式)。

对于包含 8 个 site 的 FMO 系统,总共引入了高达 411 个简谐振动自由度。在如此庞大的高维有限温度热浴下,本方法在单个变分传播轨道内便成功获得了系统的完整约化动力学映射通道。

2.3 关键计算数据分析

A. 激发态动力学与静态无序效应(图 2)

通过对 Choi 矩阵进行不同的初态收缩,作者展示了以下两个场景:

  • 情景一(图 2a):初始激子完全定域在 Site 1。结果显示了激子在非马尔可夫热浴下的相干相干转移。当加入 $\sigma = 100\text{ cm}^{-1}$ 的静态无序(Dashed lines)时,Site 1, 2, 3 的相干振荡显著被抑制,且在 600 fs 之后,Site 3 的占据数几乎趋向平稳。这说明静态无序极大地加速了局域化过程。
  • 情景二(图 2b):初始态被相干激光脉冲制备在最低激子本征态 $|E_1\rangle$。分析表明,在此激子表象下,能级间的布居数转移相当缓慢。特别地,静态无序对激子基态的布居动力学几乎没有影响,这揭示了激子相干性对无序的鲁棒性。

B. Choi 矩阵本征值谱演化:去相干与松弛的定量区分(图 3)

利用 Choi 矩阵的本征值谱 $\{\lambda_\alpha(t)\}$ 随时间的变化,可以直接表征系统-环境纠缠的本质:

  • 在 $t \to 0$ 时,光谱由单一的本征值 $\lambda_1 \approx 1.0$ 主导,表明此时通道是近乎幺正的、纯的量子演化;
  • 在极短时间内($< 100\text{ fs}$),最大本征值骤降,而少数几个次级本征值迅速抬升。这一快速阶段对应于系统的**纯去相干(Pure Decoherence)**过程,即系统相干性在不发生能级转移的情况下向热浴耗散;
  • 在较长的时间尺度($> 200\text{ fs}$)上,更多微小的本征值逐渐浮现,形成了一条长尾巴。这一慢速阶段代表了能量松弛与重分配(Population Relaxation),表明物理状态在激子本征态之间转移。

这种不依赖任何特定初始状态和基底选择的“全景诊断”,是传统初态动力学方法绝无可能提供的。

C. 冯·诺依曼熵与有效秩(图 4)

基于归一化的 Choi 状态 $\rho_{\text{Choi}}(t) = J_\Phi(t)/d_S$,其冯·诺依曼熵定义为:

$$S(t) = -\text{Tr}[\rho_{\text{Choi}} \log_2 \rho_{\text{Choi}}]$$

有效秩(Effective Rank)定义为 $r_{\text{eff}}(t) = 2^{S(t)}$。计算结果(图 4)表明:

  • $S(t)$ 和 $r_{\text{eff}}(t)$ 在初期急剧上升,在 $300-400\text{ fs}$ 达到峰值($r_{\text{eff}} \approx 20$),对应最大通道混合度;
  • 随后,$S(t)$ 与 $r_{\text{eff}}(t)$ 呈现轻微的平缓下降趋势,而不是单调上升趋于最大值 $2 \log_2 d_S = 6$(对应完全退极化通道 $\sim 64$)。这有力地证明了,即使在长时极限下,热浴也未完全随机化系统,而是保持了显著的物理结构相干性,系统趋向于平衡态的热非极化分布。

3.1 偏迹(Partial Trace)收缩的张量网络算法

构建 Choi 矩阵的关键计算瓶颈在于对环境自由度执行快速、数值稳定的偏迹收缩。论文的核心算法实现在 Appendix A 中给出,在张量列(TT)格式下,具体的偏迹求和步骤如下:

设整体 TFD 纯化态在 TT 格式下的核心张量为 $G^{(k)}$,物理与辅助系统位于位置 $2$ 和 $1$(即排序中的 $A$ 和 $S$),环境(浴)位于位置 $3$ 至 $N_B+2$。计算偏迹矩阵 $J_\Phi(i_A, i_S; j_A, j_S)$ 的具体步骤为:

  1. 右至左扫描(Right-to-Left Sweep):从最后一个热浴格点 $t = N_B + 2$ 开始,自右向左逐步收缩环境核心张量。 对于最右侧核心,构造边界矩阵: $$E^{(N_B+2)}_{\alpha_{N_B+1}, \alpha'_{N_B+1}} = \sum_{k_{N_B}} G^{(N_B+2)}_{\alpha_{N_B+1}, k_{N_B}, 1} \bar{G}^{(N_B+2)}_{\alpha'_{N_B+1}, k_{N_B}, 1}$$
  2. 对于其余的环境节点 $t = N_B + 1, \dots, 3$,递归计算转移矩阵: $$E^{(t)}_{\alpha_t, \alpha'_t} = \sum_{\alpha_{t+1}, \alpha'_{t+1}} E^{(t+1)}_{\alpha_{t+1}, \alpha'_{t+1}} \sum_{k_t} G^{(t)}_{\alpha_t, k_t, \alpha_{t+1}} \bar{G}^{(t)}_{\alpha'_t, k_t, \alpha'_{t+1}}$$ 该过程将所有的环境核心信息压缩进了 $2 \times 2$ 的过渡矩阵 $E^{(3)}_{\alpha_2, \alpha'_2}$。
  3. 系统与辅助核的组装: 构造主系统核 $S$(位置 $2$)的张量: $$M^{(2)}_{\alpha_1, \alpha'_1; i_S, j_S; \alpha_2, \alpha'_2} = G^{(2)}_{\alpha_1, i_S, \alpha_2} \bar{G}^{(2)}_{\alpha'_1, j_S, \alpha'_2}$$ 吸收压缩后的环境矩阵 $E^{(3)}$: $$\widetilde{M}^{(2)}_{\alpha_1, \alpha'_1; i_S, j_S} = \sum_{\alpha_2, \alpha'_2} M^{(2)}_{\alpha_1, \alpha'_1; i_S, j_S; \alpha_2, \alpha'_2} E^{(3)}_{\alpha_2, \alpha'_2}$$
  4. 辅助格点(Ancilla)收缩: 构造辅助核 $A$(位置 $1$)的张量: $$M^{(1)}_{i_A, j_A; \alpha_1, \alpha'_1} = G^{(1)}_{1, i_A, \alpha_1} \bar{G}^{(1)}_{1, j_A, \alpha'_1}$$ 最终两端收缩得到完整的 Choi 矩阵: $$J_\Phi(i_A, i_S; j_A, j_S) = \sum_{\alpha_1, \alpha'_1} M^{(1)}_{i_A, j_A; \alpha_1, \alpha'_1} \widetilde{M}^{(2)}_{\alpha_1, \alpha'_1; i_S, j_S}$$

3.2 复现指南与开源链接

为了方便科研界复现这一工作,作者在 GitHub 平台及 Zenodo 数据库发布了完整的计算代码及 FMO 模型参数。

  • 数据及代码托管地址(Zenodo)DOI: 10.5281/zenodo.19855477
  • 主要依赖软件包
    • TT-Toolbox / TensorToolbox (基于 Matlab 或 Python-numpy 实现的张量网络核心库)
    • 基于 One-site TDVP 算法的波函数传播求解器
    • 利用自研的 低秩插值分解(ID) 工具箱对光谱密度进行离散化(代码及数据均已在 SI 及 Zenodo 库中提供)。

3.3 运行复现步骤(以 Python 环境为例)

  1. 环境配置: 安装支持 Tensor-Train 的底层库,例如 scikit-tt 或使用自带基于 Numpy/Scipy 编写的核心收缩算子。
  2. 光谱密度构建: 调用数据包中提供的 FMO_site_specific_SDs.mat,利用离散化工具生成物理及 Ttilde 模式下的参数列表 $g_{knm}(\beta)$ 和 $\omega_k$。
  3. TDVP 演化: 初始化极大纠缠态 $|\Omega\rangle$ 与真空热浴态 $|0\rangle$。运行单位置 TDVP 传播,设定积分步长 $\Delta t = 0.25\text{ fs}$,最大张量秩为 $220$,设置截断误差限制 $\epsilon = 10^{-6}$,总传播时长 $800\text{ fs}$。
  4. Choi 矩阵提取与后处理: 在每个时间步长调用 trace_bath_tt 函数,运行 Appendix A 的自右向左收缩程序,输出并存储 $64 \times 64$ 的 $J_\Phi(t)$。利用得到的 Choi 矩阵对任意初态进行重构,并计算冯·诺依曼熵。

4. 关键引用文献,以及你对这项工作局限性的评论

4.1 关键引用文献

  1. Redfield 1965 [1] (经典马尔可夫近似理论的基石)
  2. Tanimura & Kubo 1989 [4] (HEOM 理论的开山之作)
  3. Choi 1975 / Jamiołkowski 1972 [17, 18] (Choi 同构的数学提出)
  4. Strachan, Purkayastha & Clark 2024 [19] (将 Choi 矩阵与张量网络用于费米子浴的开创性尝试)
  5. Cerrillo & Cao 2014 [15] (转移张量法 TTM 的理论提出)
  6. Takahashi & Umezawa 1996 [31] (热场动力学 TFD 纯化方法的奠基性文献)
  7. Lubich, Oseledets & Vandereycken 2015 [49] (时变变分原理 TDVP KSL 积分方法)

4.2 本工作局限性深度剖析

尽管本工作在有限温度约化动力学映射的计算效率上取得了重大突破,但作为技术作者,我认为该方法在实际推广和应用中仍存在一些无法回避的科学与技术局限性:

  1. 系统维度的指数瓶颈($d_S$ 灾难): 该方法需要将系统与一个 Ancilla 系统 $A$ 复合,这使得复合系统空间维度增加为 $d_S^2$。在张量链中,系统与辅助系统节点的物理基底维度变为了 $d_S$。对于 FMO 等 8 site 体系($d_S = 8$),$d_S^2 = 64$ 是完全可行且轻松的。然而,若要研究更大的光合捕光超级复合物、多激子相干或多带凝聚态系统(例如 $d_S = 64$),系统节点处的物理维度将增至 $4096$,这会导致系统核与浴核之间的 TT 边界秩(Rank)急剧膨胀,进而引发严重的张量收缩维数灾难。因此,该方法无法直接外推到超大型量子主系统。
  2. 迹守恒(TP)不完全性与近似截断的内在矛盾: 正如 1.5 节所述,变分 TDVP 本质上是在割裂的切空间切丛上逼近幺正演化,模守恒并不能在变分投影后自动诱导迹守恒。尽管可以通过大幅度缩减积分步长(例如从 $1\text{ fs}$ 压缩到 $0.25\text{ fs}$)使 TP 残差变小,但这显著拉长了实际计算时间。如果对于超强耦合系统,需要使用更大的 TT 秩去抑制 TP 误差,这将在一定程度上抵消张量网络原本带来的高速度优势。
  3. 静态无序高斯相干态表示的单色局限: 利用零频模式的相干态来等效静态无序的高斯抽样,是一个极为精妙的数学技巧,但该技巧仅对严格的高斯分布或特定形式分布有效。在复杂的真实分子晶体或多晶薄膜中,由于微观局域环境的多样性,静态无序分布可能呈现出非高斯(双峰、长尾等)形态。在这种非平凡分布下,零频模式的初始制备将变得极其复杂乃至不可行,迫使研究者重新回到传统的独立采样抽样轨道上去。
  4. 长时动力学外推中马尔可夫性假设的暗礁: 在第 5 节的推导中,利用转移张量法(TTM)进行长时外推,前提是必须选择一个合适的内存截断时间 $M_{\text{cut}}$。若热浴关联函数存在极长周期的振荡衰减(如强极化子耦合或极低温度下),在 $T_c$(本研究中为 $800\text{ fs}$)之内记忆核并未完全衰减。强行进行短窗口截断长时外推,无法物理地保证外推映射的完全正性(Completely Positive, CP)和迹守恒。论文中虽然验证了 FMO 在 300 K 时可稳定外推至 3 ps 且能保持 CP,但对于超低温(如 77 K 或 4 K)和强耦合,这一外推稳定性仍有待严苛检验。

5. 补充理论推导:从 Choi 映射到记忆核与有效 Pauli 动力学

为了展示该框架在化学物理后处理和动力学简化分析中的强力应用,本节补充提供从 Choi 动力学映射数据出发,推导 Nakajima-Zwanzig(NZ)非局部记忆核以及构建低维 Pauli 速率主方程的完整数学细节。

5.1 Nakajima-Zwanzig 非本地记忆核的离散时间递归推导

在开放系统理论中,非局部 Nakajima-Zwanzig 广义主方程写为:

$$\dot{\rho}_S(t) = \int_{0}^{t} K(t - s) \rho_S(s) \, ds$$

其中 $K(t)$ 即为包含所有非马尔可夫记忆效应的广义记忆核。在映射层面上,这等价于如下 Volterra 型积分方程:

$$\dot{\Phi}(t) = \int_{0}^{t} K(t - s) \Phi(s) \, ds, \quad \Phi(0) = \mathcal{I}$$

将连续时间轴离散化为均匀格点 $t_n = n \Delta t$。若对时间导数使用向前差分公式 $\dot{\Phi}(t_n) \approx \frac{\Phi_{n+1} - \Phi_n}{\Delta t}$,对卷积项使用左矩形积分规则,则上述积分方程可以离散化为:

$$\frac{\Phi_{n+1} - \Phi_n}{\Delta t} = \Delta t \sum_{m=0}^{n-1} K_{n-m} \Phi_m, \quad n \ge 1$$

我们由此可导出记忆核 $K_n$ 的离散递归计算公式:

$$K_n = \frac{\Phi_{n+1} - \Phi_n}{\Delta t^2} - \sum_{j=1}^{n-1} K_j \Phi_{n-j}, \quad n \ge 1$$

利用在上一节基于 Choi 算符重构的 $\Phi_n$ 通道系列,我们可以非常稳定地计算出任意时间的记忆核元素 $K_{n,ab}$(对应论文中图 5 展示的实部与虚部演化数据)。这表明非马尔可夫效应在大约 $300\text{ fs}$ 之后快速消散,为利用有限记忆长度进行动力学长时外推提供了铁证。

5.2 转移张量法(TTM)长时动力学演化与外推

一旦记忆核在 $t \le M_{\text{cut}}$ 时间段内衰减到零,对于任意更长的时间 $t_n > M_{\text{cut}}$,其演化可以直接通过转移张量截断格式(Transfer Tensor Method, TTM)进行显式递推:

$$\Phi_n = \sum_{m=1}^{M_{\text{cut}}} T_m \Phi_{n-m}, \quad \rho_n = \sum_{m=1}^{M_{\text{cut}}} T_m \rho_{n-m}$$

转移张量 $T_m$ 与记忆核 $K_m$ 存在一一对应关系,其在 $800\text{ fs}$ 内计算完毕后便可封存,作为动力学外推引擎(如图 6 中外推到 3 ps 甚至更长的结果,虚线部分与直接模拟高度吻合),极大地节省了模拟极长纳秒动力学的时间。

5.3 投影 Pauli 经典速率方程的建立与最佳有效速率矩阵 $W^\star$

在很多化学动力学应用中,研究人员仅关心各 site 布居数(Population)的演化规律,而不关心相干交叉项的细节。此时,可定义一个投影算符 $\mathcal{P}$,将密度矩阵映射为其对角部分(即布居数向量 $\mathbf{p}(t)$):

$$\mathcal{P}\rho = \sum_{m=1}^{d_S} |m\rangle\langle m| \rho |m\rangle\langle m|$$

由约化映射 $\Phi(t)$ 的对角部分,可构建“布居数通道映射” $M(t)$,满足:

$$\mathbf{p}(t) = M(t) \mathbf{p}(0), \quad M(t)_{i,j} = \Phi(t)_{ii,jj}$$

为了在长时极限(去相干完成、系统处于经典跳跃态)下建立形如 $\dot{\mathbf{p}}(t) = W \mathbf{p}(t)$ 的 Pauli 经典主方程,我们需要确定一个时间无关的常数有效速率矩阵 $W^\star$。作者引入了巧妙的变分最小二乘法,在预设的时间窗 $[t_{n_0}, t_{n_1}]$ 内最小化变分残差:

$$W^\star = \arg\min_{W \in \mathcal{P}} \sum_{n=n_0}^{n_1} \omega_n \left\| \dot{M}(t_n) - W M(t_n) \right\|_F^2$$

约束条件为 $W^\star$ 必须满足概率守恒及非负跃迁速率要求:

$$W_{ij} \ge 0 \quad (i \ne j), \quad W_{jj} = -\sum_{i \ne j} W_{ij}$$

这一拟合策略从最基础的量子非马尔可夫动力学映射出发,提供了一条严丝合缝、完全自洽地推导复杂体系宏观马尔可夫速率常数的通道,架起了量子相干动力学与经典化学反应动力学之间坚实的桥梁。