来源论文: https://arxiv.org/abs/2607.05269v1 生成时间: Jul 07, 2026 06:10

线性矩阵乘积切空间中的激发谱计算与秩断层扫描技术深度解析

0. 执行摘要

在现代强关联量子多体物理与量子化学计算中,精确捕获低能激发态是一项极具挑战性的任务。传统的密度矩阵重整化群(DMRG)方法在寻找一维系统的基态方面取得了举世瞩目的成功,但对于激发态的系统性构建,其效率和数学严谨性往往受限。近年来,基于矩阵乘积态(MPS)切空间(Tangent Space)的方法为这一问题提供了极其优雅且高效的变分解决方案。

本文针对 Schmidt 和 Carusotto 的最新研究成果进行深度学术解析。该工作系统性地阐明了具有开边界条件(OBC)的有限、非均匀一维系统在线性 MPS 切空间下的变分构造。更重要的是,作者引入了一种革命性的诊断工具——粒子分辨切空间秩断层扫描(Particle-Resolved Schmidt-Rank Tomography, PRSR)。该方法通过将 MPS 的全局虚拟键度规(Bond Dimension)精细化分解为特定粒子数扇区的子块秩分布,定量回答了变分切空间在不同物理扇区中的“表达能力缺陷”(Parametric Deficiency)。

通过对典型的玻色-哈巴德(Bose-Hubbard, BH)模型进行基准测试,该方法成功复现了强关联极限下的低能荷电激发(Particle/Hole Excitations),并精准捕捉了莫特绝缘体(Mott-Insulator)到超流(Superfluid)相变在有限尺寸系统中的先驱信号(荷电能隙的闭合)。本文将从代数几何、微分几何、变分动力学以及数值计算四个维度,对该方案进行全方位的深度拆解。


1. 核心科学问题、理论基础与方法细节

1.1 核心科学问题

在一维强关联多体系统中,如何高效且保真地计算非均匀系统在开边界条件下的低能激发谱?进一步地,如何从数学上定量度量和理解一个给定的基态 MPS 所生成的切空间,在表达特定物理对称性扇区激发时的精确度上限与物理表达缺陷? 传统的切空间方法常假设系统具有平移对称性(利用热力学极限下的统一切空间 ansatz),而对于边界效应显著、含有杂质或外加势场的非均匀有限系统,其切空间的代数结构和几何约束要复杂得多。本工作正是为了攻克这一理论与实用诊断工具的空白。

1.2 理论基础:代数几何视角下的 MPS 簇(MPS Variety)

从现代代数几何的视角来看,固定键度(Bond Dimension)向量 $\mathbf{D} = (D_0, D_1, \dots, D_N)$ 的矩阵乘积态不再仅仅是一个张量网络缩并的表示,而是定义在投影希尔伯特空间 $\mathbb{P}(\mathcal{H})$ 中的一个代数簇(Algebraic Variety),记为 $\mathcal{V}_{\mathbf{D},\mathbf{d}}$。这里 $\mathcal{H} = \bigotimes_{i=1}^N \mathcal{H}_i$ 为总希尔伯特空间,其局域物理维度为 $\mathbf{d} = (d_1, \dots, d_N)$。

对于开边界条件(OBC),我们强制约束边界键度 $D_0 = D_N = 1$。仿射 MPS 缩并映射 $\Phi_{\mathbf{D},\mathbf{d}}$ 将参数空间(由一系列局域张量 $\{A^i\}$ 构成)映射到多体态空间:

$$\Phi_{\mathbf{D},\mathbf{d}} : \mathcal{A}_{\mathbf{D},\mathbf{d}} \to \mathcal{H}$$

其坐标映射写为:

$$\psi = \sum_{j_1=1}^{d_1} \dots \sum_{j_N=1}^{d_N} (A^1_{j_1} A^2_{j_2} \dots A^N_{j_N}) e_{j_1} \otimes \dots \otimes e_{j_N}$$

其中 $A^i_{j_i} \in \mathbb{C}^{D_{i-1} \times D_i}$。投影 MPS 簇 $\mathcal{V}_{\mathbf{D},\mathbf{d}}$ 即为该映射图像在 Zariski 拓扑下的闭包(Zariski Closure):

$$\mathcal{V}_{\mathbf{D},\mathbf{d}} := \overline{\text{Im}(\mathbb{P}\Phi_{\mathbf{D},\mathbf{d}})}^{\text{Zar}} \subseteq \mathbb{P}(\mathcal{H})$$

根据代数几何定理,该簇可以等价地通过张量“扁平化”(Flattening)矩阵的秩约束来刻画。定义第 $k$ 个剪切位置(Cut)的扁平化算符 $T^{(k)}$:

$$T^{(k)} : \bigotimes_{i=k+1}^N \mathbb{C}^{d_i} \to \bigotimes_{i=1}^k \mathbb{C}^{d_i}$$

则 $\mathcal{V}_{\mathbf{D},\mathbf{d}}$ 内任意多体态满足:

$$\text{rank}(T^{(k)}) \le D_k, \quad \forall k = 1, \dots, N-1$$

为了消除规范自由度(Gauge Freedom)并便于数值计算,我们采用左正则规范参数化(Left-Canonical Parametrization)。引入复 Stiefel 流形:

$$\text{St}(D_i, D_{i-1}d_i) := \{ M^i \in \mathbb{C}^{D_{i-1}d_i \times D_i} \mid (M^i)^* M^i = \text{id}_{D_i} \}$$

这里局域张量被重组为矩阵:

$$M^i = \begin{bmatrix} M^1 \\ \vdots \\ M^{d_i} \end{bmatrix}$$

其满足左正则条件 $\sum_{j_i=1}^{d_i} (M^j_{j_i})^* M^j_{j_i} = \text{id}_{D_i}$。在此规范下,左正则参数空间定义为 $\mathcal{P}_{\mathbf{D},\mathbf{d}} := \prod_{i=1}^{N-1} \text{St}(D_i, D_{i-1}d_i) \times \mathbb{P}(\mathbb{C}^{D_{N-1} \times d_N})$。

1.3 规范群、垂直与水平切空间(Vertical and Horizontal Tangent Spaces)

由于参数化映射 $\Phi_{\mathbf{D},\mathbf{d}}$ 是非单射的,其内部蕴含了规范变换群的群作用。规范群定义为:

$$\mathcal{G}_{MPS} := \prod_{i=1}^{N-1} U(D_i)$$

群元 $(G_1, \dots, G_{N-1})$ 作用在左正则参数空间上,保持物理多体态不变。因此,参数空间的切空间 $T_p\mathcal{P}_{\mathbf{D},\mathbf{d}}$ 自然地分裂为两部分:

  1. 垂直切空间(Vertical Tangent Space) $\text{ker}(D_p\Phi_{\mathbf{D},\mathbf{d}})$:对应于规范轨道的切方向,这些方向纯粹由规范旋转产生,不改变物理态。其复维度为 $\sum_{i=1}^{N-1} D_i^2$。
  2. 水平切空间(Horizontal Tangent Space) $\text{ker}(D_p\Phi_{\mathbf{D},\mathbf{d}})^\perp$:与规范轨道正交的方向,代表了物理上真正的变分扰动方向。这是我们用于构建变分激发的物理切空间。

命题 7(水平切向量的正交性条件): 一个切向量 $(\delta M, \delta C) \in T_p\mathcal{P}_{\mathbf{D},\mathbf{d}}$ 属于水平切空间,当且仅当以下厄米性矩阵条件满足:

$$Y_i := \sum_{j_i=1}^{d_i} (\delta M^i_{j_i})^* M^i_{j_i} - \sum_{j_{i+1}=1}^{d_{i+1}} M^{i+1}_{j_{i+1}} (\delta M^{i+1}_{j_{i+1}})^* \quad \text{是厄米矩阵 (Hermitian)}$$

对于边界 $i=N-1$:

$$Y_{N-1} := \sum_{j_{N-1}=1}^{d_{N-1}} (\delta M^{N-1}_{j_{N-1}})^* M^{N-1}_{j_{N-1}} - \sum_{j_N=1}^{d_N} C_{j_N} (\delta C_{j_N})^* \quad \text{是厄米矩阵 (Hermitian)}$$

该定理的证明(详见论文附录)巧妙地利用了 Frobenius 度规与微分几何中的流形正交投影,为消除多体动力学中的冗余度提供了坚实的代数根基。

1.4 时变变分原理(TDVP)与谱重建算法细节

基于时变变分原理(TDVP),在流形 $\mathcal{V}^{\equiv}_{\mathbf{D},\mathbf{d}}$(代表满秩 MPS 嵌入流形)上的动力学演化可以通过将精确的薛定谔方程投影到该流形的物理切空间上来实现。Dirac-Frenkel 变分条件写为:

$$\langle \delta \psi \mid (i\partial_t - H) \psi \rangle = 0, \quad \forall \delta \psi \in T_{[\psi]}\mathcal{V}^{\equiv}_{\mathbf{D},\mathbf{d}}$$

为了提取激发谱,我们考虑在基态 $\bar{\psi}$(通过 DMRG 预先优化得到,对应变分能量 $\omega_0$)附近的线性微扰:

$$\psi_\epsilon(t) = e^{-i\omega_0 t} [\bar{\psi} + \epsilon \delta \psi(t)]$$

设水平切空间的一组实基底为 $\{h_1, \dots, h_r\}$,其通过微分映射生成的物理多体态代表元为 $\Theta_a = D_{\bar{p}}\Phi_{\mathbf{D},\mathbf{d}}(h_a)$。我们将微扰态在复化的物理切空间中展开:

$$\delta \psi(t) = \sum_{a=1}^r \xi_a(t) \Theta_a = F \xi(t)$$

其中 $F = [\Theta_1 \dots \Theta_r]$ 构成了非正交的基底矩阵。代入 TDVP 拉格朗日量并展开到二阶项,可导出一个二次拉格朗日量,其对应的 Euler-Lagrange 方程给出广义特征值问题:

$$B x = \omega A x$$

其中重叠矩阵(Overlap Matrix)$A$ 与有效哈密顿矩阵 $B$ 分别为:

$$A_{ab} = \langle \Theta_a \mid \Theta_b \rangle, \quad B_{ab} = \langle \Theta_a \mid (H - \omega_0) \Theta_b \rangle$$

重叠矩阵的压缩(SVD Compression): 在实际希尔伯特空间中,水平切空间参数化生成的 $\Theta_a$ 向量可能线性相关。因此,在求解广义特征值方程之前,必须对基底进行正交化压缩。我们对 $F$ 进行隐式奇异值分解(或直接对厄米矩阵 $A$ 进行对角化 $A = V \Sigma V^\dagger$),仅保留奇异值大于阈值(例如 $10^{-10}$)的本征模,得到正交规范化物理切空间基底 $U = [u_1, \dots, u_s]$。广义特征值问题随即转化为标准厄米对角化问题:

$$B_{red} y = \omega y, \quad B_{red} = U^\dagger (H - \omega_0) U$$

本征值 $\omega$ 即为对应激发态的激发能。


2. 关键 Benchmark 体系、计算数据与性能分析

2.1 玻色-哈巴德(Bose-Hubbard)一维链基准测试

为了评估线性 MPS 切空间方法的精度和表达缺陷,作者选用了极具代表性的一维玻色-哈巴德模型。其哈密顿量为:

$$H = -J \sum_{i=1}^{N-1} (a^\dagger_i a_{i+1} + \text{h.c.}) + \frac{U}{2} \sum_{i=1}^N n_i(n_i - 1) - \mu \sum_{i=1}^N n_i + \eta \sum_{i=1}^N (a^\dagger_i + a_i)$$
  • 参数设置:局域玻色子最大占据数限制为 $n_i \in \{0, 1, 2\}$(即局域物理维度 $d = 3$)。固定化学势 $\mu = 0.5$,相互作用强度 $U = 1$。系统尺寸 $N = 10$,采用开边界条件(OBC)。
  • 物理对称性:当 $\eta = 0$ 时,系统具有 $U(1)$ 粒子数守恒对称性,总粒子数算符 $\hat{M} = \sum_i n_i$ 为守恒量。基态与激发态均可严格归入特定粒子数 $M$ 扇区。当 $\eta = 0.01$ 时, $U(1)$ 对称性显式破缺。

2.2 计算数据与能谱重建($N=10$, $D=8$)

图 3 展示了在全局键度限制 $\mathbf{D} = [1, 3, 8, 8, 8, 8, 8, 8, 8, 3, 1]$ 下,切空间重建谱(crosses)与完全确切对角化(Exact Diagonalization, ED, dots)的对比:

  1. 对称性守恒极限($\eta = 0$)

    • 在极小 $J$ 区域(莫特绝缘相极限,基态处于 $M_0 = 10$ 扇区),重建激发谱与 ED 极其完美地吻合。最接近零能的物理激发是位于 $M = 11$ 扇区的粒子型(Particle-like)荷电激发,以及位于 $M = 9$ 扇区的空穴型(Hole-like)荷电激发。
    • 随着跳跃振幅 $J$ 的增大,由于动能项占优,荷电能隙开始收缩。在 $J_c \simeq 0.18$ 附近,系统的基态发生交叉,从 $M_0 = 10$ 跃变到 $M_0 = 11$ 扇区。切空间方法极其敏锐且准确地捕捉到了荷电模的“软化”(Softening),这是有限尺寸系统中莫特绝缘体到超流体相变的典型先驱信号。
    • 同时,基态的一体密度矩阵最大本征值 $\lambda_0$ 在接近相变点时显著增大(图 3 中部),表明非局域相干性的增强,切空间方法在如此剧烈的物理跃变区仍保持了极高的精度。
  2. 对称性破缺极限($\eta = 0.01$)

    • 由于对称性破缺,量子态无法再按粒子数分类,能级交叉演化为“避免交叉”(Avoided Crossing)。基态平均粒子数 $\langle \hat{M} \rangle$ 随 $J$ 连续光滑变化。切空间方法(图 3 底部)在这种强对称性破缺场景下,其重建低能激发谱依然与 ED 的连续演化曲线保持了惊人的一致。

2.3 精度与键度(Bond Dimension)依赖性分析

图 4 定量刻画了粒子型和空穴型激发模式在不同最大键度 $D \in \{4, 5, 6, 7, 8, 9\}$ 下的变分误差 $|\omega - \omega_{exact}|$:

  • 在深莫特绝缘体区($J \ll 1$),由于基态几乎没有纠缠(接近直积态),即使是非常小的键度(如 $D=4$),切空间方法得到的激发能误差也能达到 $\sim 10^{-11}$ 的机器精度极限。
  • 随着 $J$ 的增大,多体纠缠迅速建立,激发能误差呈现单调上升趋势。增加键度 $D$ 能显著压低误差,例如在 $J=0.15$ 处,将 $D$ 从 $4$ 提升到 $9$ 可以将误差降低近 4 个数量级。
  • 精度跃变(Precision Jumps):在 $J_c \simeq 0.18$ 的基态交叉点处,误差曲线出现明显的非连续性阶跃(表现为突然掉落变精确)。这直接暴露了基态粒子数扇区改变引起的“切空间表达能力突变”。为了深刻理解这一不连续性,必须引入文章的核心创新——秩断层扫描技术。

3. 粒子分辨施密特秩分布(PRSR)与切空间断层扫描

这是本研究最具学术深度的理论贡献:为什么在全局键度 $D$ 保持恒定时,切空间对激发谱的重建精度会在不同区域发生剧烈突变?

3.1 粒子分辨施密特秩分布(PRSR)的数学定义

考虑一个具有 $U(1)$ 对称性的总粒子数为 $M_0$ 的基态多体态 $\bar{\psi} \in \mathcal{H}_{M_0}$。选择任意一个系统剪切位置(Cut)$\ell \in \{1, \dots, N-1\}$。希尔伯特空间在此位置可双分(Bipartite)为左块和右块:

$$\mathcal{H} = \mathcal{H}^{[1:\ell]} \otimes \mathcal{H}^{[\ell+1:N]}$$

由于总粒子数守恒,左、右子空间进一步分解为局域粒子数子空间:

$$\mathcal{H}^{[1:\ell]} = \bigoplus_m \mathcal{H}^{[1:\ell]}_m, \quad \mathcal{H}^{[\ell+1:N]} = \bigoplus_n \mathcal{H}^{[\ell+1:N]}_n$$

由于总数必须为 $M_0$,基态的纠缠谱(Schmidt Decomposition)只能包含满足 $m + n = M_0$ 的配对扇区。因此,跨越剪切线 $\ell$ 的基态状态矢可精确表示为:

$$\bar{\psi}_\ell = \sum_{m=m_{min}^{(\ell)}}^{m_{max}^{(\ell)}} \bar{\psi}_\ell^{(m)}$$

其中,第 $m$ 粒子扇区对应的子块系数矩阵为 $\Psi_\ell^{(m)}$。我们定义该子块矩阵的秩为粒子分辨施密特秩(PRSR)

$$r_\ell(m) := \text{rank}(\Psi_\ell^{(m)})$$

显然,全局普通键度 $D_\ell$ 与 PRSR 的关系为:

$$D_\ell = \sum_{m} r_\ell(m)$$

核心启示:即使两个不同的量子态在全局具有完全相同的化学键度 $D_\ell$,它们内部的 PRSR 空间分布 $r_\ell(m)$ 也可能完全不同,这直接导致了由它们生成的物理切空间在不同物理激发扇区中的表达本领天差地别。

3.2 局域张量的选域规则(Selection Rules)

PRSR 的分布直接约束了 MPS 局域矩阵 $M^j_n$ 的块结构。利用 PRSR 展开,局域张量变为分块矩阵:

$$M^j_n(m, m') : \mathbb{C}^{r_j(m')} \to \mathbb{C}^{r_{j-1}(m)}$$

由于局域格点注入了 $n$ 个物理粒子,虚拟指数处的累计粒子数必然满足粒子数选择定则(Particle-Flow Selection Rule)

$$M^j_n(m, m') = 0, \quad \text{unless } m' = m + n$$

这一选域规则同样组织了切空间的一阶变分扰动 $\delta M^j_n(m, m')$。与基态张量必须死死锁定在 $m' = m + n$ 轨道不同,通用的物理切空间扰动可以探查所有满足 $r_{j-1}(m) \ne 0$ 且 $r_j(m') \ne 0$ 的虚拟块组合。当 $m' \ne m + n$ 时,这些扰动就生成了带有粒子/空穴型荷电性质的物理激发。

3.3 参数化表达能力缺陷(Parametric Deficiency)

我们定义投影到特定物理粒子数扇区 $M$ 的水平切空间为 $T^{proj}_M := P_M T$。该空间的维度可通过对角化重叠矩阵直接得到:

$$\text{dim}(T^{proj}_M) = \text{rank}(P_M U) =: \rho_M$$

我们期待的“目标激发维度”为 $d^{target}_M = \text{dim}(\mathcal{H}_M)$ (若排除了基态自身,则为 $\text{dim}(\mathcal{H}_M)-1$)。据此定义切空间在 $M$ 扇区的参数缺陷度

$$\delta^{par}_M := d^{target}_M - \rho_M$$

如果 $\delta^{par}_M = 0$,说明切空间完美覆盖了该物理扇区的所有变分自由度,谱重建将是绝对精确的(达到机器精度);若 $\delta^{par}_M > 0$,意味着切空间在该扇区存在“表达盲区”,重建激发谱必然存在近似变分误差。

3.4 $N=4, D=5$ 微观模型下的断层扫描实例剖析

为了以极高的学术透彻度展示断层扫描的魔力,我们对 $N=4$、全局键度 $D=5$ 的 Bose-Hubbard 链进行剖析(对应图 5 结果):

随着 $J$ 的变化,我们通过 PRSR 诊断识别出了四个截然不同的代数区域。其中心剪切线 $\ell=2, 3$ 上的基态 PRSR 分布、以及切空间在不同粒子数 $M$ 下的实际物理维度 $\rho_M$(本征值 $\omega$ 重建谱)对比如下:

区域 (Region)剪切 $\ell=2$ 的 PRSR $r_2(m)$剪切 $\ell=3$ 的 PRSR $r_3(m)$切空间在各扇区的实际维度 $(\rho_M)_{M=0}^8$目标维度 $d^{target}_M$物理效应分析
橙色区 (Orange)$(0, 2, 1, 2, 0, 0)$$(0, 0, 1, 1, 1, 0)$$(0, 4, 6, 16, 12, 16, 6, 4, 0)$$M=4$ 扇区: $19$$M=3,5$ 扇区: $16$对 $M=3,5$ 扇区,由于 $\rho_3=\rho_5=16=d^{target}$,参数缺陷为0,激发谱完美重建。对 $M=4$ 扇区,由于 $\rho_4=12 < 19$,参数缺陷为7,导致 $M=4$ 的高能变分激发(图 5 中部)无法被切空间捕获,出现明显的重建误差。
绿色区 (Green)$(0, 1, 2, 2, 0, 0)$$(0, 0, 1, 1, 1, 0)$$(0, 3, 8, 14, 14, 13, 8, 3, 0)$同上基态 PRSR 的局域移动,直接引起了切空间支撑维度的重组。$\rho_4$ 提升到 14,但仍然不完备。
红色区 (Red)$(0, 1, 2, 1, 1, 0)$$(0, 0, 1, 1, 1, 0)$$(1, 3, 8, 13, 14, 13, 8, 3, 1)$同上开启了全局边界粒子通道($M=0,8$ 的 $\rho$ 变为 $1$)。
蓝色区 (Blue)$(0, 1, 2, 1, 1, 0)$$(0, 0, 0, 1, 1, 1)$$(0, 2, 6, 12, 16, 13, 10, 4, 1)$同上粒子分布进一步向高占有态漂移,重组了切空间的整体拓扑分布。

这一极其精美的断层扫描结果,完美解释了为什么在图 5 中,$M=3$ 和 $M=5$ 扇区的低能与高能激发均能被计算蓝十字完全重合地复现,而 $M=4$ 扇区的高能激发分支在变分上直接“丢失”了(因为存在高达 7 个维度的参数缺失)。这绝非计算程序的数值不稳定,而是由于基态的 PRSR 代数结构对切空间维度产生的内禀代数几何约束


4. 代码实现细节、算法复现与计算复杂度

为了方便科研人员在主流张量网络库(如 ITensors.jl, TensorKit.jlMPSKit.jl)中高保真地复现该算法,本节给出核心计算流程的伪代码及算法架构设计。

4.1 算法流程伪代码设计

#= 
  线性 MPS 切空间激发谱重建核心算法
  输入: 
    Hamiltonian MPO: H_MPO
    Max Bond Dimension: D_max
    局域物理维度: d
  输出:
    激发能谱: ω_spectra
=#

using LinearAlgebra
using TensorOperations

function compute_excitation_spectra(H_MPO, D_max, d, N)
    # 步骤 1: 运行标准的 DMRG 算法获得基态 MPS (左正则化形式)
    # gs_mps 包含局域 Stiefel 矩阵 M[1..N-1] 和 右端张量 C
    gs_mps = run_dmrg(H_MPO, D_max)
    omega_0 = compute_expectation_value(gs_mps, H_MPO)
    
    # 步骤 2: 构建水平切空间的投影矩阵与实基底 {h_a}
    # 利用命题 7 的厄米条件构建约束矩阵 R 的零空间 (Nullspace)
    R_matrix = construct_tangent_constraint_matrix(gs_mps)
    # 提取实 nullspace 作为实水平参数空间基底
    H_basis = nullspace(R_matrix) 
    r = size(H_basis, 2) # 切空间实维度
    
    # 步骤 3: 隐式构建重叠矩阵 A 和有效哈密顿矩阵 B
    # 为避免直接在大希尔伯特空间中构建 q 维的 F 矩阵,利用 MPO 缩并技术直接计算元素
    A = zeros(ComplexF64, r, r)
    B = zeros(ComplexF64, r, r)
    
    println("开始构建投影算符... 总参数通道数: ", r)
    for a in 1:r
        theta_a = apply_differential_map(gs_mps, H_basis[:, a])
        for b in a:r
            theta_b = apply_differential_map(gs_mps, H_basis[:, b])
            
            # 直接在线性空间缩并收缩 (利用 MPS 快速内积)
            A[a, b] = dot(theta_a, theta_b)
            B[a, b] = dot(theta_a, H_MPO - omega_0 * id, theta_b)
            
            # 厄米共轭填充
            A[b, a] = conj(A[a, b])
            B[b, a] = conj(B[a, b])
        end
    end
    
    # 步骤 4: 执行 SVD 压缩,过滤重叠矩阵的奇异值以消除冗余
    # 对重叠矩阵进行本征对角化以确定物理独立方向
    valA, vecA = eigen(Hermitian(A))
    keep_indices = findall(x -> x > 1e-10, valA)
    s = length(keep_indices)
    println("经过冗余剪裁后的物理独立切空间维度: ", s)
    
    # 构建变换矩阵 U_transform
    U_proj = vecA[:, keep_indices] * diagm(1.0 ./ sqrt.(valA[keep_indices]))
    
    # 步骤 5: 将 B 投影到正交化的物理子空间,求解标准特征值问题
    B_reduced = Hermitian(U_proj' * B * U_proj)
    ω_spectra, y_vectors = eigen(B_reduced)
    
    return ω_spectra
end

4.2 算法复杂度与计算瓶颈分析

直接在大希尔伯特空间中构建完整的密集矩阵 $F \in \mathbb{C}^{q \times r}$(其中 $q = d^N$ 是物理维度,呈现指数爆炸)是不可行的。本算法的精妙之处在于所有的计算和收缩都在 MPS/MPO 局部链上进行

  • 物理参数数量 $P$

    $$P = \sum_{j=1}^N d_j D_{j-1} D_j$$

    切空间实维度 $r \le 2P$。对于均匀链,其复杂度为 $\mathcal{O}(N d D^2)$。

  • 收缩一对物理切空间矢量的计算复杂度: 如果哈密顿量表示为虚拟键度为 $w$ 的 MPO,收缩两个切向量对 $\langle \Theta_a \mid H' \mid \Theta_b \rangle$ 的占优转移动作(Transfer Matrix Action)复杂度为:

    $$\mathcal{O}(N w d_{\max}^2 D_{\max}^4)$$
  • 全局矩阵组装总复杂度: 由于共有 $\mathcal{O}(r^2) \sim \mathcal{O}(N^2 d^2 D^4)$ 个矩阵元素对需要计算,构建完整的 $A$ 和 $B$ 矩阵的计算总复杂度标度为:

    $$\mathcal{O}(N^3 w d^4 D^8)$$

    虽然将原本针对全希尔伯特空间的指数级 $\mathcal{O}(d^2N)$ 复杂度降为了关于键度的多项式级,但其对键度高阶项 $\mathcal{O}(D^8)$ 的强依赖性是该算法在数值应用上的核心瓶颈。这也正是为什么在此类计算中,必须引入极为高效的局部三阶张量(Three-Center)收缩算法以压低比例常数。

4.3 开源工具推荐与集成方案

  • ITensor (C++/Julia):目前国际上最流行且极度易用的强关联张量网络计算框架。可以通过其内置的 MPSMPO 数据结构以及自动指标匹配功能(Out-of-the-box Index Matching),极易编写上述物理切空间投影函数。
  • MPSKit.jl:一个专门针对 Julia 语言开发、面向切空间物理(包括 TDVP 和 VUMPS 算法)的高端开源库。其内部已高度集成了零空间构建、规范化投影以及准粒子激发谱(Quasiparticle Excitations)的构建工具,是复现本工作断层扫描的最佳起点。

5. 关键引用文献与局限性评述

5.1 关键引用文献

  1. Schollwöck (2011) [arXiv:1008.3477]:阐明了 MPS 与 DMRG 基础理论的圣经级综述,为本论文规范群 $G_{MPS}$ 的引入奠定了物理基础。
  2. Haegeman et al. (2011) [Phys. Rev. Lett. 107, 070601]:首次将时变变分原理(TDVP)引入一维晶格 MPS 切空间的开创性工作。
  3. Haegeman, Mariën, Osborne, and Verstraete (2014) [J. Math. Phys. 55, 021902]:从微分几何(Metric, Parallel Transport and Curvature)视角系统化定义 MPS 簇几何属性的奠基石。

5.2 局限性评述与学术批判

尽管该工作在有限系统非均匀切空间的代数理解上取得了突破,但在量子化学和材料物理的实际落地应用中,仍面临以下关键局限:

  1. 高阶多粒子激发的失效(Rank Frustration): 线性切空间本质上是一阶变分扰动(First-Order Linear Perturbation)。正如本工作第 3 节断层扫描所揭示的,当目标激发的粒子数偏离基态粒子数较远(例如涉及双玻色子激发、双激子激发等)时,切空间的维度 $\rho_M$ 会受到基态 PRSR 纠缠秩的极其严厉的惩罚(即 $\delta^{par}_M$ 极其巨大)。这导致线性切空间方法在处理高度非线性、多体激发过程时完全无能为力。为了克服这一缺陷,必须开发高阶变分方法(如文中展望的“双切空间变分架构”),但这将导致计算复杂度进一步呈现灾难性的飙升。
  2. 向高维物理流形的拓展极难(PEPS Obstruction): 将该理论拓展到二维晶格系统的投影纠缠配对态(PEPS)是凝聚态理论学家的夙愿。然而,PEPS 参数空间的规范结构(Gauge Structure)不是简单紧致的 Lie 群,而是具有极高代数复杂度的规范非局域闭包,其对应的“水平切空间”零空间投影算符在数学上极难严格求解,计算复杂度会发生质的改变。
  3. 大键度下的 $\mathcal{O}(D^8)$ 标度壁垒: 对于强纠缠物理系统,基态 DMRG 的键度通常需要达到 $D \sim 1000$ 甚至更高。此时 $\mathcal{O}(D^8)$ 的谱重建标度(即使在 $N=10$ 的极小系统上)将遭遇难以逾越的算力壁垒。在量子化学的高精度能谱计算中,如何引入近似收缩和局部有效近似以降低该指数,是该方法走向实用的关键步阀。

6. 补充内容:从微观矩阵实例深度理解 PRSR 机制

为了给致力于自主实现代码的科研人员提供最直观的数学脚手架,我们在此彻底拆解论文附录 C 中给出的一个 4 格点、粒子数 $M_0 = 4$ 的微观一维系统实例,看粒子分辨施密特秩如何控制具体的矩阵结构。

6.1 物理系统的状态表示

设系统尺寸 $N = 4$,局域物理基底为 $e_0, e_1, e_2$(每个格点最多占据 2 个玻色子),总粒子数 $M_0 = 4$。基态多体波函数选为:

$$\bar{\psi} = |00\rangle_L |22\rangle_R + |01\rangle_L |12\rangle_R + 2|10\rangle_L |21\rangle_R + |02\rangle_L |02\rangle_R + |11\rangle_L |11\rangle_R + |20\rangle_L |20\rangle_R$$$$+ |12\rangle_L |01\rangle_R + |12\rangle_L |10\rangle_R + |21\rangle_L |01\rangle_R - |21\rangle_L |10\rangle_R + |22\rangle_L |00\rangle_R$$

我们选择在系统中心位置 $\ell = 2$ 进行剪切,左右块基底根据粒子数进行整理。例如:

  • 左块在粒子数 $m=1$ 扇区的基底为:$\mathcal{B}_L(1) = \{ |01\rangle_L, |10\rangle_L \}$
  • 右块在粒子数 $n=3$ 扇区的基底为:$\mathcal{B}_R(3) = \{ |12\rangle_R, |21\rangle_R \}$

6.2 施密特块矩阵的分步拆解与秩计算

我们将基态按照左块包含的粒子数 $m \in \{0, 1, 2, 3, 4\}$ 进行分解,构建各个扇区的系数矩阵 $\Psi_2^{(m)}$,这些矩阵即为式 (9) 扁平化算符的直接映射:

  1. $m=0$ 块(左块 0 粒子,右块 4 粒子): 左基底 $\{|00\rangle_L\}$,右基底 $\{|22\rangle_R\}$。对应项为 $|00\rangle_L |22\rangle_R$。系数矩阵为:

    $$\Psi_2^{(0)} = \begin{bmatrix} 1 \end{bmatrix} \implies r_2(0) = \text{rank}(\Psi_2^{(0)}) = 1$$
  2. $m=1$ 块(左块 1 粒子,右块 3 粒子): 左基底 $\{|01\rangle_L, |10\rangle_L\}$,右基底 $\{|12\rangle_R, |21\rangle_R\}$。提取相关项:$|01\rangle_L |12\rangle_R + 2|10\rangle_L |21\rangle_R$。系数矩阵为:

    $$\Psi_2^{(1)} = \begin{bmatrix} 1 & 0 \\ 0 & 2 \end{bmatrix} \implies r_2(1) = \text{rank}(\Psi_2^{(1)}) = 2$$
  3. $m=2$ 块(左块 2 粒子,右块 2 粒子): 左基底 $\{|02\rangle_L, |11\rangle_L, |20\rangle_L\}$,右基底 $\{|02\rangle_R, |11\rangle_R, |20\rangle_R\}$。相关项为 $|02\rangle_L |02\rangle_R + |11\rangle_L |11\rangle_R + |20\rangle_L |20\rangle_R$。系数矩阵为:

    $$\Psi_2^{(2)} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \implies r_2(2) = \text{rank}(\Psi_2^{(2)}) = 3$$
  4. $m=3$ 块(左块 3 粒子,右块 1 粒子): 左基底 $\{|12\rangle_L, |21\rangle_L\}$,右基底 $\{|01\rangle_R, |10\rangle_R\}$。相关项为 $|12\rangle_L |01\rangle_R + |12\rangle_L |10\rangle_R + |21\rangle_L |01\rangle_R - |21\rangle_L |10\rangle_R$。系数矩阵为:

    $$\Psi_2^{(3)} = \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \implies r_2(3) = \text{rank}(\Psi_2^{(3)}) = 2$$
  5. $m=4$ 块(左块 4 粒子,右块 0 粒子): 左基底 $\{|22\rangle_L\}$,右基底 $\{|00\rangle_R\}$。相关项为 $|22\rangle_L |00\rangle_R$。系数矩阵为:

    $$\Psi_2^{(4)} = \begin{bmatrix} 1 \end{bmatrix} \implies r_2(4) = \text{rank}(\Psi_2^{(4)}) = 1$$

6.3 结论统一

将这五个粒子扇区的信息组装,我们可以得到跨越中心对称切线的 PRSR 空间特征分布

$$(r_2(m))_{m=0}^4 = (1, 2, 3, 2, 1)$$

对应的全局普通 Schmidt 键度为:

$$D_2 = \sum_{m=0}^4 r_2(m) = 1 + 2 + 3 + 2 + 1 = 9$$

这个具体而微的代数过程,极其直观地揭示了粒子守恒对称性如何深刻地约束并格式化多体状态在虚拟指数空间的分布拓扑。也正是因为这种分块结构的存在,基态才像模盘一样,将切空间的激发维度严格圈定在特定的“容许轨道”内。当系统参数 $J$ 驱动 PRSR 谱发生跃变时,切空间在特定物理扇区的有效变分基底数量(即实际维数 $\rho_M$)就会出现剧烈的非连续阶跃,这就最终对激发能谱的精度曲线(图 4)产生了决定性的断层级调制。本工作揭示的这一底层物理和代数几何联系,将极大地丰富未来一维乃至更高维强关联体系变分激发谱算法的设计与诊断方法论。