来源论文: https://arxiv.org/abs/2606.04936v1 生成时间: Jun 05, 2026 00:52
0. 执行摘要
在强关联电子系统（如高温超导体、莫特绝缘体和多自由度量子化学体系）的计算中，多体微扰论（MBPT）发挥着至关重要的作用。其中，Parquet方程（镶嵌方程）作为一种完全自洽的双粒子自洽图表方法，能够同时且无偏地处理电荷、自旋、轨道等多个物理通道的涨落及其相互纠缠。然而，长期以来，自洽Parquet迭代在进入中强耦合区域时，极易遭遇“误导性收敛”（Misleading Convergence），即算法表面上收敛，却收敛到一个完全非物理（Unphysical）的亚稳态，或者干脆彻底发散。
维也纳工业大学（TU Wien）的 Herbert Eßl 和 Stefan Rohshap 等人在其最新工作中，系统性地攻克了这一难题。他们从动力学系统的分岔理论出发，系统分析了戴森-施温格方程与贝特-萨尔皮特方程联立迭代的雅可比矩阵（Jacobian Matrix）谱特征，首次给出判定物理不动点失稳的显式代数判据。令人惊讶的是，研究证明误导性收敛不仅发生于双粒子不可约顶点发散区，在无顶点发散的参数区间同样会由于本征值穿过零点而失稳。基于此发现，他们提出了一种控制本征方向符号反转的稳定化算法，并构造了在强耦合极限下天然稳定的**“强耦合迭代”方案**。在零维模型（Zero-Point Model）和哈巴德原子（Hubbard Atom）这两大严苛的非微扰基准体系上，该算法成功跨越多条发散线，收敛到了物理分支。这一成果为将Parquet自洽计算、动力学顶点近似（D$\Gamma$A）推广至真实高维强关联格点计算扫清了最大的算法障碍。
1. 核心科学问题，理论基础，技术难点与方法细节
1.1 核心科学问题：何为“误导性收敛”？
在多体物理自洽计算中，我们通常将一系列非线性方程写成不动点迭代形式：
$$\\Psi^{(n+1)} = \\mathcal{F}[\\Psi^{(n)}]$$其中，$\\Psi$ 包含了格林函数、自能以及各类双粒子顶点函数。由于物理方程的高度非线性（例如贝特-萨尔皮特方程中顶点与格林函数的二次耦合），算符 $\\mathcal{F}$ 在参数空间（如库仑排斥能 $U$、化学势 $\\mu$）中往往存在多个不动点分支。其中仅有一个对应于物理格林函数（即物理不动点 $\\Psi^*_{\\text{phys}}$），其余均为由于代数多值性引入的非物理不动点 $\\Psi^*_{\\text{unphys}}$。
当相互作用力 $U$ 较小时，常规微扰论初值可以保证迭代收敛到 $\\Psi^*_{\\text{phys}}$。但随着 $U$ 的增大，系统发生非微扰相变或发生强关联效应，物理不动点发生分岔，其对应的不动点迭代算符的雅可比矩阵最大本征值模长超过 $1$，导致物理不动点变为排斥子（Repeller），而非物理不动点可能恰好变为吸引子（Attractor）。此时，常规迭代会“无可挽回”地滑向非物理分支，且数值误差极其顺滑地收敛至收敛阈值以下（即误导性收敛），使研究者得出完全错误的物理结论。
1.2 理论基础：Parquet方程组的数学架构
Parquet架构的基础是将完全双粒子顶点 $F$ 按照通道可约性（Reducibility）进行严格分解。对于 $SU(2)$ 对称系统，通常定义密度通道（$d$）、磁性通道（$m$）、单态配对通道（$s$）和三态配对通道（$t$）。核心方程组由以下四部分耦合而成：
贝特-萨尔皮特方程（Bethe-Salpeter Equations, BSE）： 将双粒子可约顶点 $\\Phi_r$ 与双粒子不可约顶点 $\\Gamma_r$ 联系起来（以密度通道 $d$ 为例）：
$$\\Phi_d^{\\nu\\nu'\\omega} = \\frac{1}{\\beta^2} \\sum_{\\nu_1\\nu_2} \\Gamma_d^{\\nu\\nu_1\\omega} \\chi_{0, \\text{ph}}^{\\nu_1\\nu_2\\omega} F_d^{\\nu_2\\nu'\\omega}$$其中 $\\chi_{0, \\text{ph}}^{\\nu\\nu'\\omega} = -\\beta G(\\nu) G(\\nu+\\omega) \\delta_{\\nu\\nu'}$ 为裸磁化率（泡图）。
Parquet严格分解方程： 完全不可约顶点 $\\Lambda$、各通道可约顶点 $\\Phi_r$ 与完全不可约顶点 $\\Gamma_r$ 之间的几何代数重组：
$$\\Gamma_d^{\\nu\\nu'\\omega} = \\Lambda_d^{\\nu\\nu'\\omega} - \\frac{1}{2} \\Phi_d^{\\nu(\\nu+\\omega)(\\nu'-\\nu)} - \\frac{3}{2} \\Phi_m^{\\nu(\\nu+\\omega)(\\nu'-\\nu)} + \\frac{1}{2} \\Phi_s^{\\nu\\nu'(-\\omega-\\nu-\\nu')} + \\frac{3}{2} \\Phi_t^{\\nu\\nu'(-\\omega-\\nu-\\nu')}$$施温格-戴森方程（Schwinger-Dyson Equation, SDE）： 将自能 $\\Sigma(\\nu)$ 与双粒子完全顶点 $F$ 严格联系：
$$\\Sigma(\\nu) = \\frac{U}{\\beta} \\sum_{\\nu'} G(\\nu') - \\frac{U}{2\\beta^2} \\sum_{\\nu'\\omega} (F_d^{\\nu\\nu'\\omega} - F_m^{\\nu\\nu'\\omega}) G(\\nu') G(\\nu'+\\omega) G(\\nu+\\omega)$$戴森方程（Dyson Equation）：
$$G(\\nu) = \\frac{1}{G_0^{-1}(\\nu) - \\Sigma(\\nu)}$$
在标准迭代流程中，输入裸格林函数 $G_0$、相互作用 $U$ 以及完全不可约顶点 $\\Lambda$（通常取为裸顶点 $U$ 或由DMFT给出的局域完全不可约顶点），迭代更新状态向量 $\\Psi = (\\Phi_d, \\Phi_m, \\Phi_s, \\Phi_t, G)$。为了抑制震荡，实践中常采用阻尼迭代（Damped Iteration）：
$$\\Psi^{(n+1)} = p \\mathcal{F}[\\Psi^{(n)}] + (1-p)\\Psi^{(n)}, \\quad p \\in (0, 1]$$1.3 技术难点：为什么常规阻尼迭代必然会失稳？
我们考虑不动点 $\\Psi^*$ 附近的一阶微扰线性响应。定义 $\\Pi = 1 - \\frac{\\delta \\mathcal{F}}{\\delta \\Psi}$，则阻尼迭代的雅可比矩阵可表示为：
$$J = \\frac{\\delta \\Psi^{(n+1)}}{\\delta \\Psi^{(n)}} = 1 - p \\Pi$$设 $\\lambda$ 为矩阵 $\\Pi$ 的本征值，则雅可比矩阵 $J$ 的相应本征值为 $1 - p\\lambda$。不动点局部渐进稳定的充要条件是 $J$ 的所有本征值模长均小于 1，即 $|1 - p\\lambda| < 1$。在复平面上，这一条件构成了以 $1/p$ 为半径、以 $1/p$ 为圆心的圆盘。随着相互作用 $U$ 增大，本征值 $\\lambda$ 演化，其轨迹可能通过以下三种方式越出稳定圆盘：
- 类型 I：$\\lambda$ 演化出单极点（Odd Pole），导致 $\\text{Re}(\\lambda)$ 瞬间切换正负号并趋于无穷。这通常直接对应于双粒子不可约顶点发散（Vertex Divergences）。
- 类型 II：实本征值 $\\lambda$ 连续穿过 0（即 $\\text{Re}(\\lambda) < 0$）。
- 类型 III：复共轭本征值对的实部穿过 0（即 $\\text{Re}(\\lambda) < 0$ 且 $\\text{Im}(\\lambda) \\neq 0$）。
对于类型 II 和类型 III，由于 $\\text{Re}(\\lambda) < 0$，对应的雅可比本征值 $1 - p\\lambda$ 的实部大于 1，其模长在任何正阻尼系数 $p > 0$ 下都绝对大于 1。这意味着通过一味减小阻尼系数 $p$ 无法稳定该不动点。这是阻尼迭代在强关联区失效的根本物理与数学原因。
1.4 方法细节：雅可比本征方向选择性反转算法（稳定化策略）
为了在 $\\text{Re}(\\lambda) < 0$ 时强制系统收敛至物理不动点，本工作巧妙利用了相似变换，将不稳定的本征方向对应的阻尼系数强制修正为负阻尼（即反转更新方向），而保持稳定本征方向的阻尼系数为正。具体步骤如下：
- 在迭代的第 $n$ 步，将高维状态向量 $\\Psi^{(n)}$ 展平。计算并显式构造解析雅可比矩阵 $\\Pi$。
- 对 $\\Pi$ 进行对角化，求出相似变换矩阵 $U_{\\text{eig}}$ 使得 $\\Pi = U_{\\text{eig}} \\Lambda_{\\Pi} U_{\\text{eig}}^{-1}$，其中 $\\Lambda_{\\Pi}$ 为本征值对角矩阵。
- 构造对角稳定化矩阵 $\\mathcal{D}$，其对角元按以下规则定义（其中 $p$ 为根据常规稳定条件计算的自适应步长）： $$\\mathcal{D}^{\\alpha\\alpha} = \\begin{cases} p & \\text{if } \\text{Re}(\\lambda_{\\alpha}) > 0 \\ -p & \\text{if } \\text{Re}(\\lambda_{\\alpha}) \\le 0 \\end{cases}$$ 或者使用平滑的连续映射版本： $$\\mathcal{D}^{\\alpha\\alpha} = \\begin{cases} \\delta_{\\alpha\\alpha} & \\text{if } c\\frac{2\\text{Re}(\\lambda_{\\alpha})}{|\\lambda_{\\alpha}|^2} > 1 \\ -\\delta_{\\alpha\\alpha} & \\text{if } c\\frac{2\\text{Re}(\\lambda_{\\alpha})}{|\\lambda_{\\alpha}|^2} < -1 \\ c\\frac{2\\text{Re}(\\lambda_{\\alpha})}{|\\lambda_{\\alpha}|^2} \\delta_{\\alpha\\alpha} & \\text{else} \\end{cases}$$
- 将阻尼系数对角阵变换回原始表象，构造稳定化对角矩阵（Stabilization Matrix）： $$\\mathcal{P} = U_{\\text{eig}} \\mathcal{D} U_{\\text{eig}}^{-1}$$
- 执行矩阵-向量乘积形式的更新： $$\\Psi^{(n+1)} = \\mathcal{P} \\cdot \\mathcal{F}[\\Psi^{(n)}] + (1 - \\mathcal{P}) \\cdot \\Psi^{(n)}$$
通过该动力学修正，修正后迭代系统的雅可比矩阵 $\\tilde{J} = 1 - \\mathcal{P}\\Pi$ 所有的本征值在物理不动点处的模长被强制拉回到 1 以内，从而将排斥子变为局部局部吸引子。
1.5 强耦合迭代（Strong-Coupling Iteration）方案
作为稳定化策略的备选替代，论文提出了另一项极具颠覆性的方案：迭代完全顶点 $F$ 而非可约顶点 $\\Phi$。其代数出发点在于反转贝特-萨尔皮特方程。在 $F$ 的表象下，不可约顶点可以表示为：
$$\\Gamma_r[F] = \\frac{F_r}{1 \\mp \\chi_{0,r} F_r}$$此时，BSE被严格积算，而迭代映射直接作用于 $F$。通过分析此方案下 $U=0$ 处的雅可比矩阵 $\\Pi^F$，发现其本征值恰好为 $\\lambda^F = (-1/2, 1/2, 1, 1, 1)$。由于存在实部为负的本征值 $-1/2$，该方案在弱耦合区天然不稳定，但随着 $U \\to \\infty$，其本征值全部变为正数，因而在强耦合区天然稳定。这与常规Parquet迭代的稳定性表现出了完美的互补性。
2. 关键 Benchmark 体系、计算所得数据与性能数据
为了系统检验稳定化自洽方法的有效性，论文精细研究了两个具有非微扰解析解（或高精度数值解）的基准模型。
2.1 零维模型（Zero-Point Model）
零维模型是一种不具频率和动量依赖性的简化模型，其电子仅有自旋自由度。该模型能够以最纯粹的形式展现卢廷格-沃德泛函（LWF）的多值性分支和顶点发散特征。
2.1.1 稳定性相图
论文在参数空间 $U - \\text{Re}(\\delta\\mu)$ 上绘制了自洽Parquet迭代的精确稳定性相图（见论文 Fig. 4）：
- 绿色区域（$U < 4.5$ 左右）：常规无阻尼迭代（$p=1$）稳定。所有雅可比本征值处于单位圆内。
- 黄色区域：常规阻尼迭代可稳定（属于类型 II，但可通过减小阻尼因子 $p$ 强制收敛）。
- 红色区域：常规阻尼迭代发生本质失稳，必须使用符号反转稳定化算法。
- 顶点发散点：在粒子-空穴对称点（$\\text{Re}(\\delta\\mu) = 0$）处，不稳定的边界（红绿交界）正好对应于密度通道双粒子不可约顶点的发散线（$U \\approx 9.87$, 标为黑点）。然而，偏离粒子-空穴对称性后（$\\text{Re}(\\delta\\mu) > 0$），失稳边界迅速下移至相互作用更低处（$U \\approx 4.5$），而此处根本没有顶点发散。这有力证明了顶点发散并非迭代失稳的唯一源头。
2.1.2 迭代收敛结果对比
沿着 $\\text{Re}(\\delta\\mu) = 0$ 路径（见论文 Fig. 5）：
- 当 $U < 9.87$ 时，常规迭代与稳定化迭代均收敛于物理格林函数和顶点分支。
- 越过 $U = 9.87$（第一发散线）后，常规迭代（黄色圆圈）立即偏离物理分支（灰色实线），发生了“误导性收敛”，收敛到了无物理物理意义的分支上。
- 采用本文的稳定化算法（蓝色十字），即使在 $U = 20$ 的极强相互作用下（此时已跨越了多个发散和分岔区），依然能以机器精度重构真实的非微扰物理分支（如 $\\text{Re}(\\Phi_d)$, $\\text{Im}(G)$ 等）。
2.2 哈巴德原子（Hubbard Atom, HA）
哈巴德原子是格点哈巴德模型在跃迁矩阵元 $t \\to 0$ 时的强耦合极限。它引入了无穷维松原频率（Matsubara Frequency）依赖性，因而是检验频域截断和高维数值计算稳定性的试金石。
2.2.1 频率盒子尺寸 $N_f, N_b$ 的定标与不稳定性阶数（Multiplicity）
当系统跨越顶点发散线时，不稳定的雅可比本征值数量（具有负实部的本征值个数）会发生不连续的阶跃。论文揭示了不稳定性阶数与费米松原频率截断大小 $N_f$ 的精细代数关系（见论文 Fig. 8 和 Fig. 21）：
- 在半填充（$\\delta\\mu = 0$）处，随着相互作用越过第一顶点发散线 $\\beta U_1 = 2\\pi/\\sqrt{3} \\approx 3.628$：
- 如果只考虑密度-密度子通道（Density-Density Sub-Jacobian），不稳定的本征值个数精准等于 $N_f$。
- 越过第二发散线后，不稳定本征值个数阶跃至 $2N_f$。
- 这背后的物理原因是，在顶点发散处，双粒子不可约顶点 $\\Gamma_d$ 的单发散极点通过雅可比矩阵中的克罗内克 $\\delta_{\\nu, \\tilde{\\nu}}$ 因子（由裸磁化率和格林函数代入）在频域发生了卷积扩散，导致整整 $N_f$ 个费米频率本征方向同时失稳。
- 玻色频率盒子 $N_b$ 的决定性角色：研究表明，为了使稳定化迭代能够跨越第 $n$ 条发散线，玻色子频率截断 $N_b$ 必须大于等于特定阈值（例如跨越第一条发散线要求 $N_b \\ge 5$，第二条要求 $N_b \\ge 7$，见 Fig. 26）。若玻色盒子过小，发散极点携带的频域关联信息会在通道投影变换中丢失，导致稳定化算法失效。这一重要发现为后续格点Parquet计算的频域网格划分提供了明确的量化标准。
2.2.2 算法效率：稳定化迭代 vs 拟牛顿法（Chord Method）
论文在 Fig. 12 中对比了不同自洽方案收敛到物理不动点所需的迭代步数：
- 常规阻尼方法（Uniform Damping）：在 $U > 3.6$ 之后步数剧增并最终收敛到错误分支（因此在 Fig. 12 中无数据）。
- 拟牛顿 Chord 法：在可收敛区间内，Chord 法展现出超线性的收敛速度，其迭代步数（约为 $10^2$ 次）显著低于稳定化迭代（约 $10^3 - 10^5$ 次）。
- 局限性：由于 Chord 法将所有不动点（不论物理与否）均处理为局部稳定，其收敛分支完全依赖于初值的选取。当物理不动点与非物理不动点距离极近时，Chord 法极易收敛到错误分支。而稳定化迭代通过雅可比分析，选择性地仅使物理不动点稳定，从而对初值展现出极强的鲁棒性。
3. 代码实现细节与复现指南
3.1 核心算法流图与状态向量平坦化
在编写代码时，首要难点在于将多维频域张量 $\\Phi_r^{\\nu\\nu'\\omega}$（维数大小为 $4 \\times N_f \\times N_f \\times N_b$）和单粒子格林函数 $G(\\nu)$（大小为 $N_f$）融合成一个一维平坦状态向量 $\\Psi$：
$$\\Psi = [\\vec{\\Phi}_d, \\vec{\\Phi}_m, \\vec{\\Phi}_s, \\vec{\\Phi}_t, \\vec{G}]$$平坦化后，雅可比矩阵 $\\Pi$ 的总维度为 $N_{\\text{state}} \\times N_{\\text{state}}$，其中 $N_{\\text{state}} = 4 N_f^2 N_b + N_f$。
+-----------------------+
| 输入 G0, U, Lambda |
+-----------+-----------+
|
v
+-----------------------+
| 建立初始物理猜测 Psi^0 |
+-----------+-----------+
|
v
+-----------------------+ <---+
| 求解标准 Parquet 映射 | |
| F_new = F[Psi^n] | |
+-----------+-----------+ |
| |
v |
+-----------------------+ |
| 计算解析雅可比矩阵 Pi | |
+-----------+-----------+ |
| |
v |
+-----------------------+ |
| diagonalization: | |
| Pi = U * diag(L) * U^-| |
+-----------+-----------+ |
| |
v |
+-----------------------+ |
| 构建选择性阻尼算符 D | |
| Re(L)<=0 ? D=-p : D=p | |
+-----------+-----------+ |
| |
v |
+-----------------------+ |
| P = U * D * U^-1 | |
+-----------+-----------+ |
| |
v |
+-----------------------+ |
| 执行选择性符号更新: | |
| Psi^n+1 = P*F_new + | |
| (1-P)*Psi^n | |
+-----------+-----------+ |
| |
v |
/\\ \\ |
/ \\ 收敛 (err<1e-6)? \\ |
\\ / -------------------> 结束
\\/ 否
3.2 雅可比子块解析求导矩阵元（密度通道示例）
要实现解析雅可比，必须根据附录 B 给出的解析导数公式进行矩阵元装配。以下给出密度通道映射对密度通道可约顶点的导数矩阵元装配逻辑：
$$\\frac{\\delta f_d^{\\nu\\nu'\\omega}}{\\delta \\Phi_d^{\\tilde{\\nu}\\tilde{\\nu}'\\tilde{\\omega}}} = -\\frac{1}{\\beta^2} \\left[ -\\frac{1}{2} \\delta_{\\nu,\\tilde{\\nu}} \\delta_{\\nu+\\omega,\\tilde{\\nu}'} \\chi_{0,\\text{ph}}^{(\\nu+\\tilde{\\omega})\\nu_2\\omega} (\\Gamma_d^{\\nu_2\\nu'\\omega} + \\Phi_d^{\\nu_2\\nu'\\omega}) + \\Gamma_d^{\\nu\\nu_1\\omega} \\chi_{0,\\text{ph}}^{\\nu_1\\nu\\tilde{\\omega}} \\left( \\delta_{\\nu',\\tilde{\\nu}'}\\delta_{\\omega,\\tilde{\\omega}} - \\frac{1}{2}\\delta_{\\tilde{\\nu}+\\omega,\\tilde{\\nu}'}\\delta_{\\nu'-\\nu,\\tilde{\\omega}} \\right) \\right]$$在编码时，为避免高维循环导致的极低效率，必须利用 Julia 或 Python (NumPy) 的张量缩约函数（如 opt_einsum 或 @tensor 宏）进行重组。
3.3 Julia 核心复现参考代码
以下是基于 Julia 语言，利用相似变换动态调整雅可比矩阵更新方向的稳定化求解核心代码段：
using LinearAlgebra
"""
stabilized_step!(Psi, F_map, G0, U, p; c=0.1)
执行一步自适应符号反转的稳定化更新步。
- `Psi`: 当前的一维展开状态向量
- `F_map`: 计算 Parquet 映射的函数, 返回与 Psi 维度相同的向量
- `Jacobian_generator`: 计算解析雅可比矩阵 Pi = I - dF/dPsi 的函数
"""
function stabilized_step!(Psi::Vector{ComplexF64}, F_map::Function, Jacobian_generator::Function, p::Float64; c=0.1)
N = length(Psi)
# 1. 计算当前的映射值
Psi_new = F_map(Psi)
# 2. 生成雅可比矩阵 Pi (Pi = I - dF/dPsi)
# 注意：这里的 Pi 对应于论文中的 \$\Pi = 1 - \\delta\\mathcal{F}/\\delta\\Psi\$
Pi = Jacobian_generator(Psi)
# 3. 本征值分解
# 由于 Pi 通常是非对称复矩阵，使用 eigen 进行广义对角化
eig_decomp = eigen(Pi)
vals = eig_decomp.values
U_eig = eig_decomp.vectors
U_inv = inv(U_eig)
# 4. 根据本征值实部，自适应确定对角阻尼矩阵 D
D = zeros(ComplexF64, N, N)
for i in 1:N
lambda = vals[i]
# 计算局部稳定性圆盘条件：基于公式 (13) 和 (27)
if real(lambda) > 0
# 稳定/常规阻尼区
damping_factor = min(1.0, c * 2.0 * real(lambda) / (abs2(lambda)))
D[i, i] = damping_factor
else
# 不稳定本征方向：强制进行负阻尼反转更新
damping_factor = min(1.0, c * 2.0 * abs(real(lambda)) / (abs2(lambda)))
D[i, i] = -damping_factor
end
end
# 5. 重构物理表象下的稳定化矩阵 P = U * D * U^-1
P = U_eig * D * U_inv
# 6. 执行矩阵-向量形式的混合更新
# Psi^(n+1) = P * Psi_new + (I - P) * Psi
Psi_next = P * Psi_new + (I - P) * Psi
# 计算当前步的残差
residual = norm(Psi_next - Psi) / norm(Psi)
return Psi_next, residual
end
3.4 推荐开源软件包
目前，学术界已有复现该论文及开展相关研究的开源工具。核心推荐使用：
- ParquetIR.jl (Julia 开发，对应文献 [54])：
- Repo 链接：https://github.com/S-Badr/ParquetIR.jl
- 简介：该包采用稀疏中间表象（Intermediate Representation, IR）高效存储和求解多频 Parquet 方程，内置了高度优化的解析雅可比组装模块和自洽求解器，极适合作为复现本稳定化算法的底座。
- NLsolve.jl (Julia)：用于拟牛顿法（如 Broyden, Chord）的非线性求解器，可与上述稳定化代码进行收敛性能基准对比。
4. 关键引用文献及工作局限性微评
4.1 关键引用文献
本工作建立在近年来对多体计算发散性、多值性规律的研究基础之上，以下文献是理解本项目背景的基石：
- [1] H. Eßl, M. Reitner, E. Kozik, and A. Toschi, Phys. Rev. Lett. (2026)（doi:10.1103/zjy7-4jqd）： 评价：这是本工作的前序理论篇，首次揭示了自洽微扰论中误导性收敛的雅可比根源，本工作则是将其方案在双粒子自洽Parquet体系上的具体落地。
- [5] E. Kozik, M. Ferrero, and A. Georges, Phys. Rev. Lett. 114, 156402 (2015)（doi:10.1103/PhysRevLett.114.156402）： 评价：开创性地指出自洽粗体微扰论中 Luttinger-Ward 泛函的多值性及可能收敛到非物理分支的现象（“强关联计算的梦魇”）。
- [6] R. Rossi and F. Werner, J. Phys. A: Math. Theor. 48, 485202 (2015)（doi:10.1088/1751-8113/48/48/485202）： 评价：正式提出零维模型（ZP Model），将其作为严格分析自洽微扰论收敛性的数理沙盒。
- [25] T. Schäfer, G. Rohringer, O. Gunnarsson, S. Ciuchi, G. Sangiovanni, and A. Toschi, Phys. Rev. Lett. 110, 246405 (2013)（doi:10.1103/PhysRevLett.110.246405）： 评价：首次在多体计算中发现双粒子不可约顶点的物理发散线，指出了微扰自洽方案在强关联区的失效边界。
4.2 局限性深度微评
尽管本工作在数理逻辑和稳定性控制上取得了突破，但在走向真实材料的实用计算时，仍存在以下三大不可回避的局限：
极高的计算复杂度与内存开销（Dimensionality Bottleneck）： 稳定化算法需要显式组装和对角化雅可比矩阵 $\\Pi$。对于哈巴德原子，矩阵大小为 $(4N_f^2 N_b + N_f) \\times (4N_f^2 N_b + N_f)$。若取常规中等频域格网 $N_f=12, N_b=13$，矩阵维度已达 $\\sim 7500$。如果将其应用到具有动量 $q$ 和频域依赖性的真实格点计算（如 2D Square Lattice $16 \\times 16$ 尺寸），状态向量维度将暴增至 $10^9$ 以上，显式组装和对角化雅可比矩阵（复杂度 $O(N^3)$）将彻底变得在物理上不可行。这一“维度灾难”使得目前该稳定化方案仅能用于局域或杂质模型。
局部稳定与全局吸引域的矛盾（Local vs. Global Basin of Attraction）： 雅可比稳定化算法本质上是一种局部控制。它能够保证，一旦迭代轨迹进入物理不动点附近的微扰半圆形邻域（ Basin of Attraction），迭代就绝不会逸出。但是，它无法保证从极远初值（如任意微扰初值或极弱排斥初值）出发的全局迭代轨迹一定会跨越重重非物理不动点的重围，进入该物理邻域。实践中仍需要十分精巧的物理延拓（如对相互作用 $U$ 进行极其细致的参数扫描 continuation），这极大地增加了计算的时间和人工成本。
顶点发散点处的数学奇异性（Discontinuity at Vertex Divergence Line）： 在不稳定性类型 I（即不可约顶点 $\\Gamma$ 发生单极点发散）的位置，可约顶点 $\\Phi$ 同样会发生发散。这意味着物理参数线在穿过发散线时是代数不连续的。即使稳定化算法能在两侧稳定，也无法直接在数学发散点上直接求解，这给参数 continuation 计算在临界点附近的稳定性带来了严峻考验。
5. 补充技术细节与前沿展望
5.1 雅可比本征值的物理内涵：它代表了什么？
我们不应仅仅将雅可比矩阵 $\\Pi$ 的本征谱视为一个“纯数学工具”，实际上它与系统的物理磁化率和涨落响应有着深刻的联系。
在物理上，贝特-萨尔皮特方程的本征值直接对应于系统在特定通道中的物理响应（极化率）。当某一物理极化率 $\\chi_r$ 趋向于无穷（如磁不稳定性或超导配对不稳定性临界点）时，对应的 BSE 算符的分母趋向于零，直接导致雅可比矩阵 $\\Pi$ 出现极小本征值（接近零）。因此，迭代不稳定的发生，往往是物理涨落剧烈程度的数值映像。类型 III 的复共轭失稳，在物理上可能对应着系统内不同物理通道（如反铁磁涨落与超导配对涨落）之间的强烈竞争和动力学杂化。
5.2 前沿救赎之路：QTT（量子化张量列）与单玻色子交换（SBE）
针对上述局限性1（维度灾难），学术界目前正沿着两条前沿路线对该算法进行优化提升，可望在不久的将来实现高维格点模型的稳定自洽：
基于量子化张量列（Quantics Tensor Trains, QTT）的雅可比压缩： 文献 [23, 31] 指出，多体格林函数和顶点函数在松原频域具有极强的指数级信息冗余。QTT 技术能够将大小为 $N_f^2 N_b$ 的顶点矩阵压缩为位宽仅为对数级 $\\log(N)$ 的张量网络。如果在 QTT 压缩空间中直接定义和求算雅可比算符，则其对角化开销有望从 $O(N^3)$ 直接断崖式下跌至 $O(\\log^3 N)$，从而攻克高维矩阵求解的瓶颈。
单玻色子交换（Single-Boson Exchange, SBE）表象下的稳定化： SBE（文献 [21, 79]）将完全四点顶点 $F$ 解耦为通过动力学玻色子中介的三点费米子-玻色子顶点（Fermion-Boson Vertex, $\\gamma$）以及非相互作用的三点标量。在 SBE 表象下，状态向量的阶数直接从四点（$N_f^2 N_b$）降至三点（$N_f N_b$）。在这一表象下构建稳定化雅可比，矩阵维度将下降 2-3 个数量级，极大地拓展了该算法在真实关联材料计算中的实用空间。