来源论文: https://arxiv.org/abs/2606.28855v1 生成时间: Jul 04, 2026 10:01

非平衡电子-声子动力学的高动量分辨率模拟:热化瓶颈与声子色散效应的 QTT-NEGF 深度解析

0. 执行摘要

在凝聚态物理和超快泵浦-探测(Pump-Probe)实验中,光激发固体材料后的非平衡电子-声子(e-ph)动力学是决定材料热化、退相干及诱导非平衡相变(如瞬态超导、电荷密度波重构)的核心物理过程。传统的理论模拟面临着严峻的挑战:波函数方法(如 ED、DMRG)受限于声子希尔伯特空间的无限维度;而非平衡格林函数(NEGF)方法虽然能自然处理无限声子态,但其双时间(Two-time)传播特性带来了致命的 $O(t_{\max}^2)$ 内存瓶颈和 $O(t_{\max}^3)$ 或 $O(t_{\max}^4)$ 的计算复杂度。这导致以往的研究不得不引入“局域自能”近似(如动态平均场理论 DMFT),从而完全忽略了声子色散及动量依赖的电声相互作用反馈。

瑞士弗里堡大学的 Maksymilian Środa 与 Philipp Werner 在其最新工作中,展示了**量子张量列非平衡格林函数(Quantics Tensor Train Non-Equilibrium Green’s Function, QTT-NEGF)**框架在解决上述微观动力学物理问题中的革命性威力。作者利用 QTT 对双时间传播子进行指数级压缩,在自洽重正化 Migdal 近似下,实现了高达 $256 \times 256$ 个动量格点的自洽二维非平衡电声系统模拟。该研究对光学声子与声学声子体系进行了系统性对比,首次揭示了电声热化过程中的分层瓶颈效应(Hierarchy of relaxation bottlenecks):

  1. 光学声子体系:除了确认经典的“声子能量能窗”(Phonon-window)瓶颈外,首次揭示了一个由于电子-空穴自发湮灭导致的、被截半的“还原能窗”(Reduced window $\mathcal{W}_{1/2} = [-\Omega_r/2, \Omega_r/2]$)。该能窗内呈现持久的费米子亏缺(Deficit),而窗外则存在持久的过剩(Excess);同时揭示了由于声子与粒子-空穴连续体解耦导致的“声子热化瓶颈”(Phonon thermalization bottleneck)。
  2. 声学声子体系:尽管声学声子在 $q \to 0$ 极限下具有任意小的激发能,但由于能量与动量同时守恒的约束,声子能量能窗依然顽固存在,并具有强烈的动量依赖性。由于低动量电声耦合强度的消失($g_{\mathbf{q}} \sim \sqrt{q} \to 0$)和散射的强方向不对称性,低动量声学声子模式呈现极度缓慢的发射特征,形成了极其持久的弛豫瓶颈。

本技术博客将面向具有量子化学、凝聚态物理和计算材料学背景的研究人员,对这一前沿工作的科学物理图景、数学公式推导、数值实现细节、复现步骤及未来发展方向进行深度技术解析。


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

1.1 核心科学问题:电声热化的微观动量瓶颈与局域近似的失效

当固体系统被超快激光泵浦后,电子系统瞬间吸收能量形成非热化的高能分布。这些“热”电子主要通过发射声子将能量传递给晶格,从而实现系统的整体热化。然而,这一能量和动量交换过程绝非简单的“两温度模型”(Two-Temperature Model, TTM)所能描述:

  • 非热化分布与能窗阻滞(Phonon-window bottleneck):若电子系统激发能低于单个声子能量 $\Omega$,或者电子弛豫到费米面附近的残留能量不足以发射一个声子,同时由于泡利不相容原理的限制,电子将无法向下跃迁,导致弛豫过程发生停滞,在费米面附近的一定能量范围内形成长寿命的非平衡“阻塞区”。
  • 声子重正化反馈(Renormalization feedback):声子并不是一个恒定的热库。在真实的体系中,声子通过与电子系统的自洽耦合会被剧烈地重正化(软化、寿命缩短甚至发生电荷密度波失稳)。处理这一双向反馈(Mutual feedback)必须自洽地求解电子与重正化声子的耦合方程。
  • 动量色散(Dispersion)的重要性:实际材料中的声子具有明显的动量依赖色散关系(如声学支 $\Omega_{\mathbf{q}} \sim q$,光学支在不同动量处具有不同展宽)。局域近似(如 DMFT 自能 $\Sigma(\omega)$ 与动量无关)会彻底抹杀动量守恒限制、van Hove 奇异性导致的散射各向异性,以及低动量长波声子的发射特性。要想真正描述二维或三维真实材料中的超快动力学,必须在非平衡态下全面保留动量分辨率。

1.2 理论模型基础:Fröhlich 哈密顿量与色散模型

论文采用二维平方晶格(格点数 $L \times L$)上的半满 Fröhlich 哈密顿量进行建模,在时域上的形式为:

$$H(t) = H_e + H_{ph} + H_{e-ph}(t)$$

其中,电子动能项为紧束缚近似(只考虑最近邻跃迁 $h$):

$$H_e = \sum_{\mathbf{k}\sigma} \epsilon_{\mathbf{k}} c^{\dagger}_{\mathbf{k}\sigma} c_{\mathbf{k}\sigma}, \quad \epsilon_{\mathbf{k}} = -2h(\cos k_x + \cos k_y)$$

自由声子项($\mathbf{q}$ 为动量):

$$H_{ph} = \sum_{\mathbf{q}} \Omega_{\mathbf{q}} b^{\dagger}_{\mathbf{q}} b_{\mathbf{q}}$$

电声相互作用项在 $t \ge 0$ 时被突然开启(Quench 协议):

$$H_{e-ph}(t) = \frac{1}{\sqrt{N_{\mathbf{k}}}} \sum_{\mathbf{k}\mathbf{q}\sigma} g_{\mathbf{q}}(t) c^{\dagger}_{\mathbf{k}+\mathbf{q},\sigma} c_{\mathbf{k},\sigma} (b^{\dagger}_{-\mathbf{q}} + b_{\mathbf{q}})$$

为了对比不同的色散特征,定义了三种 bare 声子色散 $\Omega_{\mathbf{q}}$ 模型:

  1. 光学声子(Optical): $$\Omega_{\mathbf{q}} = \Omega_0$$
  2. 声学声子(Acoustic): $$\Omega_{\mathbf{q}} = \frac{\Omega_0}{\sqrt{2}} \sqrt{\sin^2 \frac{q_x}{2} + \sin^2 \frac{q_y}{2}}$$
  3. 过渡情况(Intermediate): $$\Omega_{\mathbf{q}} = x \Omega_0 + (1-x) \frac{\Omega_0}{\sqrt{2}} \sqrt{\sin^2 \frac{q_x}{2} + \sin^2 \frac{q_y}{2}}$$

电声相互作用矩阵元随动量的演变必须物理上合理。对于声学声子,在长波极限($q \to 0$)下,形变势耦合导致 $g_{\mathbf{q}} \propto \sqrt{q}$。因此设计耦合系数的形式为:

$$g_{\mathbf{q}}(t) = \sqrt{\Omega_{\mathbf{q}}} g(t)$$

在 $t < 0$ 时,耦合常数 $g(t) = 0$,系统处于温度 $\beta=10$(以紧束缚带宽 $W=8h$ 为能量基准,设 $h=1$)的无相互作用热平衡态;在 $t=0$ 时,突然淬灭(Quench)开启耦合,使 $g(t \ge 0) = g = 0.4$。系统的总能量在淬灭后守恒,最终将热化到一个温度更高的自洽平衡态。

1.3 技术难点:非平衡 Kadanoff-Baym 方程的双时间计算墙

在非平衡态格林函数(NEGF)框架下,我们需要求解定义在 Keldysh 轮廓线(Keldysh contour)上的自洽 Dyson 方程:

$$G_{\mathbf{k}}(t, t') = G_{0\mathbf{k}}(t, t') + \left[ G_{0\mathbf{k}} * \Sigma_{\mathbf{k}} * G_{\mathbf{k}} \right](t, t')$$$$D_{\mathbf{q}}(t, t') = D_{0\mathbf{q}}(t, t') + \left[ D_{0\mathbf{q}} * \Pi_{\mathbf{q}} * D_{\mathbf{q}} \right](t, t')$$

其中 $*$ 代表在轮廓线上的积分褶积:$[A * B](t, t') = \int_{\mathcal{C}} d\bar{t} A(t, \bar{t}) B(\bar{t}, t')$。在实际数值求解时,这些方程会被拆分为四个分量:Matsubara($M$)、Retarded($R$)、Left-mixing($\rceil$)和 Lesser/Greater($$),合称 Kadanoff-Baym 方程(KBE)。

致命难点:

  1. 双时间存储与非马尔可夫演化:格林函数 $G(t, t')$ 和自能 $\Sigma(t, t')$ 是双时间矩阵。如果时间步数为 $N_t$,存储一个动量点的格林函数需要 $O(N_t^2)$ 空间,计算积分褶积需要 $O(N_t^3)$ 运算。由于记忆效应,每一次时间推进都需要读取之前的全部双时间数据,传统的无损存储在 $N_t > 1000$ 时就会耗尽服务器内存。
  2. 动量网格乘数效应:引入声子色散和非局域自能意味着在布里渊区有 $L \times L$ 个独立的动量格点,每个动量格点都拥有一套双时间矩阵。若 $L=256$,动量格点数高达 $65536$。在传统 NEGF 中,这会导致内存需求乘以 $6.5 \times 10^4$,这在以往是绝对无法计算的超级灾难。

1.4 方法细节:重正化 Migdal 近似与 QTT 压缩技术

为了克服这一困难,该工作引入了**量子张量列(Quantics Tensor Train, QTT)**压缩技术,并在自洽的重正化 Migdal 近似下迭代求解 KBE。

自洽自能公式(实空间表示):

为了降低动量褶积的计算代价,通常将格林函数变换到实空间求解乘积,然后再变换回动量空间。自洽电声自能 $\Sigma$ 和声子极化自能 $\Pi$ 的非马尔可夫时空表达式为:

$$\Sigma_{ij}(t, t') = i G_{ij}(t, t') V_{ij}(t, t')$$$$\Pi_{ij}(t, t') = -2i G_{ij}(t, t') G_{ji}(t', t)$$

其中 $V_{ij}(t, t') = g_i(t) D_{ij}(t, t') g_j(t')$ 是声子介导的电子-电子有效相互作用。实空间与动量空间的格林函数通过二维快速傅里叶变换(FFT)相互投影:

$$G_{ij}(t, t') = \frac{1}{N_{\mathbf{k}}} \sum_{\mathbf{k}} e^{-i\mathbf{k} \cdot (\mathbf{r}_i - \mathbf{r}_j)} G_{\mathbf{k}}(t, t')$$

QTT 压缩的核心机理:

QTT 是基于矩阵张量化(Tensorization)及张量网络分解(TT/MPS)的高效无损/有损压缩方法。其核心思想是将连续(或极细致离散化)的时间坐标进行二进制拆分: 设双时间坐标格点 $t, t'$ 对应于二进制指数 $2^R$。我们构造一个射影映射,将原本维度为 $2^R \times 2^R$ 的双时间大矩阵 $A(t, t')$ 转化为具有 $2R$ 个虚拟指标的张量:

$$t = (t_1, t_2, \dots, t_R)_2, \quad t' = (t'_1, t'_2, \dots, t'_R)_2, \quad t_m, t'_m \in \{0, 1\}$$

通过将时间自变量按位交错排列成多维张量 $A(t_1, t'_1, t_2, t'_2, \dots, t_R, t'_R)$,并利用奇异值分解(SVD)将其压缩成张量列(Tensor Train)形式:

$$A(t, t') = M^{(1)}_{t_1, t'_1} M^{(2)}_{t_2, t'_2} \dots M^{(R)}_{t_R, t'_R}$$

每个张量核心(Tensor Core)的局部矩阵大小由键合维度(Bond Dimension)$D$ 决定。对于物理上平滑的、具有局部特征或多尺度分离(Scale Separation)的格林函数,QTT 压缩能够将存储复杂度从 $O(4^R)$ 降到惊人的 $O(R \cdot D^2)$。这意味着即使在超细的时间步长下(如 $dt \sim 10^{-7}$,使得时间网格数极大),存储需求也仅随对数 $R$ 线性增长。所有在线性的 Dyson 方程求解、乘积和褶积步骤都在 QTT 格式下不解压直接进行,从而彻底攻克了“双时间计算墙”。


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

2.1 Benchmark 系统配置

  • 晶格维度:二维平方晶格,$L \times L = 128 \times 128$ 以及最大至 $256 \times 256$ 个格点。
  • 物理参数:带宽 $W = 8$(设最近邻跳符 $h=1$),声子特征频率 $\Omega_0 = 1$。初态温度 $\beta = 10$(逆温度)。淬灭耦合强度 $g = 0.4$。
  • 最大传播时间:$t_{\max} = 200$。由于采用极细的时间网格(利用 QTT 的对数增长优势),其时间网格具有极高的时间分辨率,避免了传统大步长带来的伪影。

2.2 物理计算所得关键数据与机理分析

研究通过分析费米子(电子/空穴)密度偏差 $\Delta n^k(t) = n_{\mathbf{k}}^{e,h}(t) - n_{\mathbf{k}}^{e,h}(\infty)$ 与声子数偏差 $\Delta n^q(t) = n_{\mathbf{q}}^{ph}(t) - n_{\mathbf{q}}^{ph}(\infty)$ 随时间的发展,得出了以下极其丰富的物理图像:

                     +----------------------------------------+
                     |  激光超快泵浦 (Quench 激发电声相互作用)  |
                     +-------------------+--------------------+
                                         |
                   +---------------------+---------------------+
                   |                                           |
                   v                                           v
       +-----------------------+                   +-----------------------+
       |     光学声子体系      |                   |     声学声子体系      |
       +-----------+-----------+                   +-----------+-----------+
                   |                                           |
         +---------+---------+                       +---------+---------+
         |                   |                       |                   |
         v                   v                       v                   v
  +------------+      +------------+          +------------+      +------------+
  | 电子能窗   |      | 声子色散   |          | 色散能窗   |      | 低动量声子 |
  | 还原瓶颈   |      | 极化阻滞   |          | van Hove点 |      | 发射阻滞   |
  |  W_1/2     |      | (特定动量) |          | (动量依赖) |      | (g_q -> 0) |
  +------------+      +------------+          +------------+      +------------+

2.2.1 光学声子体系(无色散体系)的弛豫与全新发现的“截半还原能窗”

  1. 主能窗瓶颈 $\mathcal{W}$: 对于光学声子,重正化后的声子频率为 $\Omega_r \approx 0.9$。经典能窗 $\mathcal{W} = [-\Omega_r, \Omega_r]$ 限制了电子的单粒子弛豫。在窗外,高能电子通过快速连续发射声子进行雪崩式弛豫,这通常在 $t \lesssim 25$ 内完成。然而,一旦电子落入 $\mathcal{W}$ 内,由于声子发射后的末态在费米面以下已被占满(泡利阻挡),弛豫通道几乎完全被关闭,导致 $\mathcal{W}$ 边界处的非平衡电子形成长寿命积累。
  2. 截半还原能窗(Reduced Window $\mathcal{W}_{1/2}$): 利用高空间和动量分辨率,作者首次在自洽 NEGF 模拟中观察到了一个极其优雅的子能窗物理现象——还原能窗 $\mathcal{W}_{1/2} = [-\Omega_r/2, \Omega_r/2]$。在非平衡演化过程中,由于电子-空穴对称性,处于 $[\Omega_r/2, \Omega_r]$ 的高能电子会与处于 $[-\Omega_r, -\Omega_r/2]$ 的高能空穴发生相干湮灭并伴随声子发射。这一复合过程将中间区域 $\mathcal{W}_{1/2}$ 内的费米子迅速排空,最终在该子能窗内形成一个明显的电子/空穴赤字区(Deficit,蓝色阴影),而在 $\mathcal{W} \setminus \mathcal{W}_{1/2}$ 区间则形成了持久的过剩区(Excess,红色阴影)。由于这些残留电子/空穴缺乏匹配的可湮灭伙伴,这一不均匀结构在 $t=200$ 时依然极具生命力,形成了难以逾越的热化瓶颈。这一发现在先前不考虑重正化的无反馈模拟或非自洽近似中被完全掩盖了。
  3. 光学声子自身的弛豫瓶颈: 过去人们默认声子能快速吸收热量。然而自洽格林函数表明,声子数偏离值 $\Delta n^q(t)$ 在动量空间中呈现出一个奇特的、以 $\mathbf{q}=(0,0)$ 为中心的红外“金刚石形”冻结赤字区。其深层物理机理在于:声子自能极化率 $\Pi_{\mathbf{q}}$ 对应于林哈德(Lindhard)电荷响应函数。在长波极限 $\mathbf{q} \to 0$ 处,电荷响应被强力压制,使得声子与粒子-空穴连续体彻底解耦。由于缺乏局域自能重正化反馈(如图 3(c) 中 $\mathbf{q}=(0,0)$ 内的零虚部所示),这些低动量光学声子模式几乎完全保持 bare 状态,无法通过重吸收多余的非平衡电子来调节自身的数目,从而使得多余的非平衡声子过剩在较宽的动量空间中长期无法消散。

2.2.2 声学声子体系(强色散体系)的各向异性与低动量发射阻滞

  1. 动量依赖能窗的解析划分: 声学声子在 $\mathbf{q} \to 0$ 时能量趋向于 0,直觉上,任何微小的电子跃迁都可以通过发射极长波的声学声子来完成,因此瓶颈似乎应该消失。然而自洽高分辨率动量谱线给出了完全相反的答案:声学体系不仅存在瓶颈,而且该瓶颈具有极为复杂的空间动量依赖性。根据严格的能量守恒 $\omega_{\mathbf{k} \to \mathbf{K}} = \epsilon_{\mathbf{K}} + \Omega_{\mathbf{K}-\mathbf{k}}$ 限制,电子发射声子的最低能量界限被深刻地限定。对于 van Hove 奇异点 $\mathbf{K}=(\pi, 0)$,其能窗宽度恰好为 $\pm \Omega_0$;而对于对角线费米面点 $\mathbf{K}=(\pi/2, \pi/2)$,其能窗边界缩窄为 $\pm \Omega_0/\sqrt{2}$。这一多重边界线在图 2(e)-(h) 中被作者绘制的理论临界线精确地圈定出来,证明自洽动力学完全遵循微观守恒律。
  2. 布里渊区各向异性非对称热化: 由于色散的存在,van Hove 点附近的电子态密度极大,导致强烈的散射。高能电子会快速地从 $\mathbf{k}=(\pi, 0)$ 方向向对角线方向 $(\pi/2, \pi/2)$ 进行“单向式”(Directional)各向异性散射。在长期演化后,对角线区域由于低能阻挡积累了大量的超额电子,而 van Hove 点附近则留下了无法填补的空穴赤字。这种空间各向异性的不平衡无法通过 $\mathbf{q} \to 0$ 的散射进行消除。
  3. 低动量 $\mathbf{q}$ 声子瓶颈: 由于变形势耦合强度 $g_{\mathbf{q}} \sim \sqrt{q} \to 0$,声学声子在接近 $\mathbf{q}=(0,0)$ 时与电子系统的耦合能力呈线性衰减。同时,由于前述费米子不对称性对散射相空间的限制,低动量长波声学声子的发射效率低得惊人。这直接导致系统虽然整体能量已经极高,但靠近布里渊区中心的低 $\mathbf{q}$ 模式声子数却无法增长,在极其漫长的时间内仍处于相对“极冷”的状态(表现为图 5(g) 中 $\mathbf{q}=0$ 附近的深蓝色赤字圆盘)。这种低动量声学声子发射的极度滞后,构成了声学色散体系最顽固的非平衡特征。

2.3 性能数据:张量网络 QTT 的革命性表现

在如此大规模的 2D 晶格自洽 NEGF 演化中,QTT 展现出了无与伦比的压倒性计算性能:

计算物理参数 \ 方法特征传统无损双时间格林函数方法本工作使用的 QTT-NEGF 方法
最大支持空间动量点数$\sim 100 \times 100$ (局部自能无反馈)$256 \times 256$ (全局非局域反馈色散)
演化最大双时间尺寸$t_{\max} \approx 30$$t_{\max} = 200$ (物理演化深度提升 7 倍)
内存占用级别 (单变量)$O(N_k \cdot N_t^2) \approx$ 几十 GB 至数 TB$O(N_k \cdot R \cdot D^2) \approx$ 几百 MB (压缩率超 $99.9\%$)
自洽收敛控制精度易在长时间传播中发生数值崩塌不收敛块时间步进 (Block-stepping) 保证极高稳定性
QTT 最大键合维度 $D_{\max}$不适用$50 \sim 80$ (对于 $128^2$);$50$ (对于 $256^2$)
SVD 截断误差阈值 $\tau$不适用$10^{-11} \sim 10^{-10}$ (高保真度无物理失真)
运行硬件平台需巨型超算节点组,极难跨越内存墙单台双路 AMD EPYC 7742 节点 (128核), 768GB RAM

3. 代码实现细节、复现指南与开源软件包

为了方便科研人员在自己的物理模型中推广该 QTT-NEGF 方法,本节将详细梳理基于 Julia 语言的算法框架和数据流复现指南。

3.1 核心算法架构与模块流图

QTT-NEGF 算法的执行遵循以下严密的自洽迭代和空间变换步骤:

+-------------------------------------------------------------+
|                        初始状态设置                         |
|         利用 TCI (张量交叉插值) 产生无相互作用的            |
|     Matsubara 格林函数 $G_M$ 和实时间格林函数 $G_0$         |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|                  第一级动量网格模拟 ($64^2$)                 |
|  循环块时间步 $t_{\text{step}}$:                            |
|  1. 将 QTT 格式的 $G_{\mathbf{k}}$ FFT 变换到实空间 $G_{ij}$ |
|  2. QTT 乘积计算自能:$\Sigma_{ij} \propto G_{ij} D_{ij}$ |
|  3. FFT 变换回动量空间自能 $\Sigma_{\mathbf{k}}$             |
|  4. 构建并求解 Dyson 扫频算法,自洽更新 $G_{\mathbf{k}}, D_{\mathbf{q}}$ |
|  5. 收敛判定:差异小于 $\epsilon_{\text{conv}} = 10^{-3}$     |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|                  第二级动量网格拓展 ($128^2$)                |
|  1. 将 $64^2$ 上收敛的自能 $\Sigma_{\mathbf{k}}$ 傅里叶插值 |
|     外推到 $128^2$ 动量网格                                 |
|  2. 利用外推自能作为初猜,在 $128^2$ 上进行全局 contour 更新|
|  3. 自洽迭代收敛至 $5 \times 10^{-3}$                       |
+------------------------------+------------------------------+
                               |
                               v
+-------------------------------------------------------------+
|                最终高级动量网格拓展 ($256^2$)                |
|  1. 傅里叶插值自能 $\Sigma$ 推进至 $256^2$ 动量网格         |
|  2. 直接计算输出非平衡电子/声子密度分布,避免超高代价的全局轮廓更新|
+-------------------------------------------------------------+

3.2 关键 Julia 依赖包与开源 Repo 链接

该研究的代码生态完全建立在 Julia 现代科学计算栈上,核心依赖包如下:

  1. ITensors.jl
    • 功能:提供底层的张量网络、矩阵乘积态(MPS)、张量列(TT)以及奇异值分解(SVD)和截断算法支撑。它是构建 QTT 离散化算符的物理基石。
    • 链接https://github.com/ITensor/ITensors.jl
  2. Tensor4all 库
    • 功能:由论文合作团队开发的专用张量网络算法库,高度优化了针对两时间积分方程的扫频更新(Sweeping solver)和张量交叉插值(TCI)。
    • 链接https://tensor4all.org (内含专门针对格林函数和 KBE 方程求解的算子映射实现)
  3. 数据和视频补充材料

3.3 实战复现指南:核心超参数控制与计算逻辑

想要成功复现该计算,必须在 Julia 脚本中精准设置以下关键控制参数:

# Julia 伪代码及超参数配置示例
using ITensors
using FFTW

# 1. 物理参数设置
const L = 128                  # 动量晶格尺寸 L x L
const beta = 10.0               # 初始逆温度
const g_coupling = 0.4          # 淬灭后的电声耦合强度
const omega_0 = 1.0             # 裸声子特征频率
const t_max = 200.0             # 最大演化时间

# 2. QTT 压缩控制超参数 (至关重要!控制计算精度与内存平衡)
const tau_cutoff = 1e-10        # SVD 截断误差,低于此值的奇异值将被丢弃
const max_bond_dim = 80         # 最大允许键合维度 D_max,防爆内存的上限屏障
const dt_step = 10.0            # 块时间步进尺寸 dt_step
const epsilon_conv = 1e-3       # 自洽 Dyson 方程迭代收敛阈值

# 3. 构造二进位时间网络
# 设格点数 N_time = 2^R
# 通过逐层引入张量指标建立自洽时间格林函数 QTT 对象
# ... (调用 ITensors 接口构造二进制链)

# 4. 核心自洽迭代循环控制
function run_self_consistent_loop(G_k, D_q)
    iter = 0
    converged = false
    while !converged
        # Step A: 快速傅里叶变换到实空间
        G_real = perform_2D_fft_to_real(G_k)
        D_real = perform_2D_fft_to_real(D_q)
        
        # Step B: 计算自能 (在 QTT 格式下不解压直接做按元素点乘)
        Sigma_real = compute_qtt_electron_self_energy(G_real, D_real, g_coupling)
        Pi_real = compute_qtt_phonon_self_energy(G_real, g_coupling)
        
        # Step C: 自能变换回动量空间
        Sigma_k = perform_2D_fft_to_momentum(Sigma_real)
        Pi_q = perform_2D_fft_to_momentum(Pi_real)
        
        # Step D: 求解 Dyson 方程更新格林函数
        G_k_new = solve_qtt_dyson_electron(Sigma_k, G_k)
        D_q_new = solve_qtt_dyson_phonon(Pi_q, D_q)
        
        # Step E: 判定收敛
        diff = calculate_qtt_difference(G_k_new, G_k) / calculate_qtt_norm(G_k)
        if diff < epsilon_conv
            converged = true
            println("Loop converged with residual: ", diff)
        end
        G_k = G_k_new
        D_q = D_q_new
        iter += 1
    end
    return G_k, D_q
end

关键复现技巧提示:

  1. 千万避免中途解压张量:QTT-NEGF 的所有算力优势均建立在格林函数处于二进制 MPS 状态。绝对不要在 FFT 或自能计算时将其转化为普通的 Array{Float64,2} 矩阵,否则内存将瞬间爆满。
  2. 采用多级外推(Multigrid extrapolation):由于直接求解 $256 \times 256$ 动量格点的自洽轮廓自能计算代价高昂。务必严格按照作者提供的分级傅里叶插值流程,在小网格($64^2$)上实现自洽后,提取自能通过三角函数插值延展到大网格,作为大网格初始值直接输出。这样可以在保持自洽反馈物理正确性的同时,节省 $95\%$ 以上的计算时间。

4. 关键引用文献与局限性评论

4.1 关键引用文献解析

本项研究依托于非平衡强关联和张量网络交叉领域的几项最重要学术奠基工作:

  • 电声弛豫与能窗阻滞的经典理论
    • [13] M. Sentef, A. F. Kemper et al., Phys. Rev. X 3, 041033 (2013):首次利用非平衡态格林函数模拟证实了超快激发下存在声子能窗引起的电子弛豫瓶颈。
    • [15] Y. Murakami, P. Werner et al., Phys. Rev. B 91, 0410128 (2015):利用重正化 Migdal 局域近似,探讨了强自洽电声重正化对瓶颈弛豫的调制作用,为本工作奠定了 Migdal 近似物理框架。
  • QTT-NEGF 核心数学方法的发展
    • [70] H. Shinaoka, M. Wallerberger, P. Werner, A. Kauch et al., Phys. Rev. X 13, 021015 (2023):首次将 Quantics Tensor Train 用于压缩非平衡态双时间物理格林函数,奠定了算法的数学可行性。
    • [61] K. Inayoshi, M. Środa, P. Werner, H. Shinaoka, arXiv:2509.15028 (2025):开发了基于因果性的分块时间步进(divide-and-conquer block-stepping)算法,彻底攻克了长时间演化下的数值不稳定性。
  • 物理机制分析基石
    • [84] G. D. Mahan, Many-Particle Physics (2000):用于推导动量依赖多体自能散射率公式(式 15)和研究声子-等离激元杂化的理论工具书。

4.2 对这项工作的局限性与前沿评论

尽管这一研究成果堪称非平衡强关联动力学领域的里程碑,但作为资深技术作者,我们仍需指出其在化学和真实材料计算应用中的几项关键物理和算法局限性:

1. Migdal 近似在极强电声耦合区(Strong-coupling regime)的失效

  • 物理局限:Migdal 近似(等价于非平衡多体微扰论的二阶自能插图)严格依赖于 Migdal 定理,即电子速度远大于声速(小无量纲参数 $m/M \ll 1$),从而忽略了所有包含电声顶点校正(Vertex Corrections)和非交叉高阶自能图。然而,在诸如富氢高温超导体(如 $\text{H}_3\text{S}$ 族)、极性半导体、过渡金属硫族化合物(TMDs)等强耦合体系中,无量纲耦合常数 $\lambda > 1.5$ 甚至更大。在此强区,顶点校正不可忽略,极易形成小极化子(Polaron)自陷域。现有的 QTT-NEGF 在 Migdal 近似下无法处理极化子凝聚和电声强关联电荷物理。

2. QTT 压缩率对“物理纠缠度”长期演化的极限挑战

  • 数学局限:QTT 能够取得极高压缩比的根本原因在于双时间函数在快/慢时间尺度上的分离。然而,当系统在强相互作用下经历长期的相干非平衡激荡,或者系统存在多频杂化、准粒子长寿命相干振荡时,双时间自变量之间的多体纠缠(Entanglement)会急剧累积。这会导致 QTT 的键合维度 $D$ 随时间 $t_{\max}$ 的延长而呈现指数级上升趋势。一旦 $D > 300$,QTT 的扫频求解器计算代价将超越传统直接求解法,从而暴露出“张量网络虚假压缩”的局限。本工作最大外推至 $t=200$,其键合维度已限制在 $D=80$,若要推进到纳秒($ns$)尺度,其算法稳定性仍待检验。

3. 与第一性原理(Ab-initio)接口的鸿沟

  • 计算局限:本研究使用的哈密顿量为高度简化的单带 Fröhlich 模型。在真实的量子化学和固体系统中,材料拥有多能带(Multi-band)、各向异性非局域跃迁,以及复杂的声子多支色散(如 3 支声学支与多支光学支的混合)。利用密度泛函微扰理论(DFPT)输出的非均匀电声矩阵元 $g_{\mathbf{k},\mathbf{q}}^{\nu, mn}$ 极度庞大且复杂,如何将这套多能带 DFPT 数据无损地映射并构建成高效的 QTT 算符矩阵乘积算子,是实现材料特异性(Material-specific)超快实验复现必须跨越的天堑。

5. 补充物理图景:微观杂化、 charge response 与 tr-ARPES 实验连接

为了给科研界同仁提供更深层的直觉,本节补充讨论自洽 QTT-NEGF 得出的物理谱线与实验可观测量的直接连接。

5.1 电荷响应与重正化声子谱的深层物理杂化

本工作最精彩的理论亮点之一是图 3 展示的重正化声子色散与电荷易受性(Charge susceptibility)的共振杂化。在强电声相互作用下,声子会诱导电荷密度起伏,声子传播子 $D_{\mathbf{q}}(\omega)$ 的自能即为粒子的极化泡 $\Pi_{\mathbf{q}}(\omega)$:

$$\Pi_{\mathbf{q}}(\omega) \approx g_{\mathbf{q}}^2 \chi_{\mathbf{q}}^{charge}(\omega)$$

在图 3(a) 与 3(b) 的对比中,我们可以清晰地观察到,在布里渊区的特定边缘,声子谱线与林哈德函数所包含的粒子-空穴连续体边界发生相交。由于自洽相互作用,声子线在这里发生分叉和弯曲,形成了凝聚态物理中典型的反交叉(Anticrossing)杂化结构,原本纯净的声子振动与电子跃迁发生电荷杂化,变成了电子-声子混合激子模式。 传统的研究由于动量网格粗糙(如仅能计算 $16 \times 16$ 动量点),无法解析出这些反交叉点在动量边界处的精细弯曲结构。本工作依靠 QTT-NEGF 提供的 $256 \times 256$ 超高分辨率动量分辨率,使得这些物理反交叉线的边缘纤毫毕现,为精确计算电荷密度波(CDW)材料中的软化自洽不稳定性提供了前所未有的显微镜。

5.2 与超快时间分辨角分辨光电子能谱(tr-ARPES)的桥梁

实验物理学家该如何观测本工作中预言的“截半还原能窗 $\mathcal{W}_{1/2}$”和“声学声子发射延迟”?答案正是时间分辨角分辨光电子能谱(Time-resolved ARPES, tr-ARPES)

tr-ARPES 直接测量的物理量是动量resolved的电子非平衡占据数:

$$I_{\text{tr-ARPES}}(\mathbf{k}, \omega; t_{\text{delay}}) \propto \text{Im} G^<_{\mathbf{k}}(\omega, t_{\text{delay}})$$

其中 $G^<_{\mathbf{k}}(\omega, t_{\text{delay}})$ 为 Lesser 格林函数的局部维格纳变换。通过将本文自洽计算出的 Lesser 成分进行局域傅里叶分析,可以直接合成出高度逼真的 tr-ARPES 理论能带谱。如图 5(a) 和 5(c) 所示:

  • 实验观测标志 A(光学声子材料):在光激发后,tr-ARPES 的光谱强度会在费米面附近展现出显著的、宽度恰为 $\pm \Omega_0/2$ 的**“光谱烧孔”(Spectral burn hole)**。在这个截半还原能窗内,光谱强度会长期低于热平衡值,而在能窗外 $[\Omega_0/2, \Omega_0]$ 区间则展现出长寿命的电子超额累积,形成双层亮带。这一非平衡指纹是判定材料中存在自洽重正化电声湮灭的终极实验证据。
  • 实验观测标志 B(声学声子材料,如富氢超导体/石墨烯):在 tr-ARPES 的时间演化图谱中,van Hove 点 $(\pi, 0)$ 的光谱衰减速度将显著快于对角线方向 $(\pi/2, \pi/2)$。由于长波声学声子的阻滞,靠近费米面内部的相干准粒子能带将展现出极其缓慢的冷却(Delayed cooling),其典型衰减寿命长达数个皮秒,比常规估算慢一个数量级以上。

本工作建立的高分辨率 QTT-NEGF 方法,不仅扫清了多维非平衡动力学计算中的“两时间存储墙”,更为建立紧密连接超快泵浦-探测实验与微观强关联理论的定量计算物理桥梁铺平了道路。