来源论文: https://arxiv.org/abs/2606.16424v1 生成时间: Jun 16, 2026 19:02

突破十亿格点极限：利用张量网络计算三维非厄米拓扑材料的实空间谱函数

0. 执行摘要

在开放量子系统、光子晶体、声学超材料以及具有增益和损耗的凝聚态物理体系中，非厄米（Non-Hermitian）系统展现出了超越传统厄米体系的丰富拓扑现象，例如非厄米皮肤效应（Non-Hermitian Skin Effect, NHSE）、异常点（Exceptional Points, EPs）以及独特的更高阶拓扑角态（Higher-Order Topological Corner Modes）。然而，对这些体系进行定量研究一直面临着巨大的计算瓶颈：随着空间维度的增加，系统的希尔伯特空间（Hilbert Space）呈指数级增长。特别是在三维（3D）有限尺寸体系中，强烈的边界效应和复杂的非厄米谱性质使得传统的精确对角化（Exact Diagonalization, ED）和稀疏矩阵求解器在格点数达到 $10^4$ 级别时便难以为继，从而无法直接揭示宏观极限（Macroscopic Limit）下的本征拓扑行为。

针对这一关键科学难题，来自芬兰阿尔托大学的 Yitao Sun、Jose L. Lado 以及瑞典查尔姆斯理工大学的 Guangze Chen 组成的研究团队，于 2026 年发表了一项具有里程碑意义的工作。他们成功开发了一种全新的张量网络（Tensor Network）理论框架。该方法首次将量化张量交叉插值（Quantics Tensor Cross Interpolation, QTCI）与非厄米核多项式方法（Non-Hermitian Kernel Polynomial Method, NHKPM）有机结合，将原本难以处理的大规模三维非厄米紧束缚哈密顿量压缩为紧凑的矩阵乘积算符（Matrix Product Operator, MPO）形式。利用该算法，研究人员在超过十亿（$10^9$）格点的三维非厄米更高阶拓扑绝缘体（HOTI）中，直接计算出了实空间局部态密度（LDOS）和谱函数。

通过调节体系的损耗强度（Loss Strength $\lambda$），该工作在宏观尺度上清晰地观测到了非厄米拓扑角态在弱损耗和强损耗机制下的演化规律，并首次揭示了强损耗下角态流形（Corner-Mode Manifold）发生的层次性分裂现象（Hierarchical Splitting）。这一突破性的成果不仅证明了张量网络算法在处理超大规模非厄米量子物态上的卓越能力，也为凝聚态物理、量子化学和材料科学中大尺寸、强无序或准晶体系的谱学模拟开辟了崭新的道路。

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

1.1 核心科学问题：三维非厄米体系的尺度壁垒

在非厄米物理学中，拓扑边界态（如一维的端点态、二维的角态、三维的角态）往往对边界的几何形状和系统尺寸高度敏感。在小尺寸晶格中，强烈的有限尺寸效应（Finite-Size Effects）会严重干扰我们对拓扑本征性质的判定，甚至可能将平庸的杂质态误判为拓扑角态。为了真正触及“宏观极限”，物理学家需要模拟包含数百万乃至数亿个格点的超大体系。然而，三维体系的格点数 $N = L^3$ 随着线性尺寸 $L$ 呈三次方增长，这意味着 $1024 \times 1024 \times 1024$ 的三维晶格拥有超过 $10^9$（十亿）个自由度。面对如此庞大的希尔伯特空间，传统的数值方法存在以下根本性制约：

  1. **精确对角化（ED）**的计算复杂度为 $\mathcal{O}(N^3)$，内存消耗为 $\mathcal{O}(N^2)$，其计算极限通常局限在 $10^4$ 左右的格点数。
  2. ** Krylov 子空间方法（如 Lanczos 或 Arnoldi 方法）**在处理非厄米矩阵时极不稳定，这是由于非厄米算符的条件数（Condition Number）可能非常庞大，极易导致数值发散。
  3. **传统核多项式方法（KPM）**虽然能够将复杂度降至 $\mathcal{O}(N)$，但它仅适用于厄米算符，因为其理论基础——切比雪夫多项式展开（Chebyshev Expansion）只有在实数区间 $[-1, 1]$ 上才收敛，无法直接应用于复平面上的非厄米谱函数计算。

1.2 理论基础：Hermitization 技巧与非厄米核多项式方法（NHKPM）

为了在保持计算线性的同时处理非厄米哈密顿量 $H$ 的复数谱，本研究采用了 Hermitization 技巧。对于一个一般的非厄米哈密顿量 $H$，其谱函数（Spectral Function）定义在复能量平面 $\omega \in \mathbb{C}$ 上：

$$f(\omega) = \langle \psi_L | \delta^2(\omega - H) | \psi_R \rangle$$

其中 $\delta^2(\omega - H)$ 是复平面上的二维狄拉克 $\delta$ 函数，$\langle \psi_L|$ 和 $|\psi_R\rangle$ 分别代表系统的左、右本征态（或特定的实空间局部探针态）。为了规避复数域上的非收敛展开，研究人员构造了一个扩大的厄米哈密顿量 $\tilde{H}(\omega)$：

$$\tilde{H}(\omega) = \begin{pmatrix} 0 & \omega I - H \\ \omega^* I - H^\dagger & 0 \end{pmatrix}$$

通过这一巧妙的映射，原本非厄米哈密顿量 $H$ 的二维谱问题被转化为了关于厄米算符 $\tilde{H}(\omega)$ 的格林函数（Green’s Function）求解问题：

$$f(\omega) = \frac{1}{\pi} \partial_{\omega^*} \langle L | (E - \tilde{H})^{-1} | R \rangle \Big|_{E=0}$$

其中，扩充基矢定义为 $|L\rangle = (0, |\psi_L\rangle)^T$ 和 $|R\rangle = (|\psi_R\rangle, 0)^T$。由于 $\tilde{H}(\omega)$ 是严格厄米的，其本征值完全为实数，因此可以安全地使用切比雪夫多项式进行展开。利用第一类切比雪夫多项式 $T_n(x)$，上式可以写为无穷级数的形式：

$$f(\omega) = \frac{2}{\pi^2} \sum_{n=1}^{\infty} (-1)^{n+1} \langle L | \partial_{\omega^*} T_{2n-1}(\tilde{H}) | R \rangle$$

在实际数值计算中，该级数通过引入 Jackson 核（Jackson Kernel）进行截断和正则化，以有效抑制切比雪夫展开带来的吉布斯振荡现象（Gibbs Oscillations）。

1.3 技术难点：多项式算符的张量网络表示与演化

尽管将问题厄米化了，但直接在十亿维度的空间中进行多项式乘法计算依然是不可能的任务。这正是张量网络（Tensor Networks）大显身手之处。该方案的核心在于：将波函数表示为矩阵乘积态（Matrix Product State, MPS），将哈密顿量及切比雪夫算符表示为矩阵乘积算符（Matrix Product Operator, MPO）。

要实现这一目标，面临两个主要的技术难点：

  1. 极大规模哈密顿量的指数压缩： 如何将十亿格点的 3D 紧束缚哈密顿量转化为占用极小内存的高效 MPO？
  2. 切比雪夫递推过程中的键维度爆炸（Bond Dimension Explosion）： 随着递推步数 $n$ 的增加，算符 $T_n(\tilde{H})$ 的算符键维度（MPO Bond Dimension）通常会呈指数增长。如果直接进行算符乘法，计算会迅速崩溃。

1.4 方法细节：量化张量交叉插值（QTCI）与 MPO 递推

为了攻克上述难点，论文提出了一套高度精密的算法流程：

1.4.1 空间格点的二进制编码与拟自旋表象

首先，算法将包含 $N = 2^L$ 个格点的实空间格点索引 $\alpha$ 编码为二进制形式：

$$\alpha = (s_1, s_2, \dots, s_L), \quad s_i \in \{0, 1\}$$

这种编码将单粒子的实空间希尔伯特空间精确映射到了一个包含 $L$ 个格点的 $\frac{1}{2}$-拟自旋（Pseudo-spin）链上。在此基底下，实空间哈密顿量 $H_{\alpha\beta}$ 可以被写成一个 $L$ 场位置张量的收缩形式，即 MPO：

$$\mathcal{H}_{(s_1, \dots, s_L), (s'_1, \dots, s'_L)} = \Gamma^{(1)}_{s_1, s'_1} \Gamma^{(2)}_{s_2, s'_2} \dots \Gamma^{(L)}_{s_L, s'_L}$$

通过利用 Quantics Tensor Cross Interpolation (QTCI) 算法，复杂的实空间函数（如调制损耗空间分布 $V(x, y, z)$）可以被极其高效地压缩进 MPO 中，其所需的虚拟键维度 $m$ 往往极小（通常在 $10 \sim 30$ 之间），从而实现了空间维度的对数级（$\mathcal{O}(\log N)$）内存占用压缩。

[实空间 3D 晶格 (10^9 格点)] 
       │ (空间坐标拆分为二进制串)
       ▼
[1D 拟自旋链 (Qubit Chain)] ──(QTCI 压缩)──► [紧凑的哈密顿量 MPO (H)]
                                                 │
                                                 ▼
[构建 Hermitized MPO (H~)] ◄──(附加一维辅助张量)──┘

1.4.2 算符层面的 Hermitization

为了在 MPO 层面直接构造 $\tilde{H}(\omega)$，算法在拟自旋链的头部（最大尺度端）引入了一个额外的“辅助格点”（Auxiliary Site），将其物理维度设为 2（对应泡利算符 $\sigma^\pm$）。通过以下张量积形式构建 $\tilde{H}$：

$$\tilde{H} = \sigma^+ \otimes (\omega I - H) + \sigma^- \otimes (\omega^* I - H^\dagger)$$

这一操作避免了显式地写出大矩阵，直接在 MPO 格式下完成了算符的厄米化，且仅将 MPO 的链长增加了 1，从而完美控制了计算复杂度。

1.4.3 切比雪夫算符导数的递推关系

非厄米谱函数不仅需要计算切比雪夫多项式本身，还需要计算其对复能量 $\omega^*$ 的导数。递推关系如下：

  • 切比雪夫多项式递推：

    $$T_0(\tilde{H}) = I, \quad T_1(\tilde{H}) = \tilde{H}$$

    $$T_{n+1}(\tilde{H}) = 2\tilde{H} T_n(\tilde{H}) - T_{n-1}(\tilde{H})$$

  • 导数算符递推：

    $$\partial_{\omega^*} T_0(\tilde{H}) = 0, \quad \partial_{\omega^*} T_1(\tilde{H}) = \partial_{\omega^*} \tilde{H} = \sigma^- \otimes I$$

    $$\partial_{\omega^*} T_{n+1}(\tilde{H}) = 2 (\partial_{\omega^*} \tilde{H}) T_n(\tilde{H}) + 2 \tilde{H} (\partial_{\omega^*} T_n(\tilde{H})) - \partial_{\omega^*} T_{n-1}(\tilde{H})$$

在计算过程中，状态矢量的递推是通过 MPS-MPO 乘法（并伴随着变分奇异值分解截断，Variational SVD Truncation）来实现的。通过引入最大键维度限制（Maximum Bond Dimension $\chi$），算法能够防止多项式迭代过程中的状态纠缠度爆炸，确保整个计算的单步复杂度严格控制在 $\mathcal{O}(L \chi^3)$ 级别。

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

2.1 物理模型：三维非厄米更高阶拓扑绝缘体（3D NH-HOTI）

为了验证该算法的精确度与可扩展性，研究团队设计了一个经典的三维 Benalcazar-Bernevig-Hughes (BBH) 模型，并引入了空间调制的虚部势能（即非厄米损耗）。系统的总哈密顿量为：

$$H = H_0 + H_V$$

其中，厄米部分 $H_0$ 描述了具有二聚化跃迁（Dimerized Hopping）和交错磁通（Staggered $\pi$-flux）的三维紧束缚模型：

$$H_0 = \sum_{\mathbf{r}, \mu=x,y,z} \left[ s_\mu(\mathbf{r}) t(\mu) c^\dagger_{\mathbf{r}+\hat{\mu}} c_{\mathbf{r}} + \text{h.c.} \right]$$

交错相位结构用于在每个小方腔（Plaquette）中引入 $\pi$ 磁通，定义为：

$$s_x(\mathbf{r}) = 1, \quad s_y(\mathbf{r}) = (-1)^x, \quad s_z(\mathbf{r}) = (-1)^{x+y}$$

跃迁振幅具有交替二聚化形式：$t(\mu) = t_1$（当坐标 $n_\mu$ 为奇数时），$t(\mu) = t_2$（当坐标 $n_\mu$ 为偶数时）。当 $t_1 < t_2$ 时，该系统处于更高阶拓扑绝缘体（HOTI）相，在三维立方的 8 个顶点处支持零维的拓扑角态（Corner Modes）。

非厄米损耗部分 $H_V$ 定义为：

$$H_V = i \lambda \sum_{x,y,z} V(x, y, z) c^\dagger_{x,y,z} c_{x,y,z}$$

其中 $\lambda$ 代表整体损耗强度，空间调制函数 $V(x,y,z)$ 沿三个方向进行余弦调制：

$$V(x, y, z) = \prod_{\gamma=x,y,z} \left[ \sqrt{2} f(a_\gamma) \cos\left(\frac{\pi}{2}\gamma + \frac{\pi}{4}\right) \right]$$

这里 $a_\gamma \equiv \gamma \pmod 4$，且 $f(a_\gamma) = +1$ (若 $a_\gamma = 0, 3$)； $f(a_\gamma) = -1$ (若 $a_\gamma = 1, 2$)。该损耗项打破了体系的厄米性，重塑了复平面上的能谱，使拓扑角态从体能带（Bulk Continuum）中spectrally分离出来，从而便于通过局部实空间谱进行观测。

2.2 小体系 Benchmark 结果（格点数 $N = 16^3 = 4096$）

在小尺寸体系下，研究人员首先将该张量网络算法（NHKPM）与精确对角化（ED）的结果进行了直接对比，以确立算法的精度。计算参数设为：$t_1 = 1$, $t_2 = 1.4$。

2.2.1 弱损耗机制（Weak-Loss Regime, $\lambda = 2$ 或 $1.4$）

  • ED 能谱： 如图 2(a) 所示，在复能量平面上，体能谱分布在两侧，而在虚轴中间（$\text{Re}(E) = 0$）附近存在一个独立的拓扑角态能区（由红圈标出）。
  • LDOS 分布： 如图 2(b) 所示，此时计算实空间局部态密度，可以发现角态在实空间中表现为高度局域化，绝大部分波函数密度集中在 3D 晶格最外层的角格点上（Outermost corner）。

2.2.2 强损耗机制（Strong-Loss Regime, $\lambda = 5$）

  • ED 能谱演化： 当损耗强度增加到 $\lambda = 5$ 时，情况发生了戏剧性的变化（图 2(c)）。原本简并或集中的角态能谱分裂成了四个在复平面上相互分离的孤立能级（In-gap Modes）。
  • 层次性局域化（Hierarchical Localization）： 分别提取这四个分裂能级的实空间 LDOS（图 2(d)-(g)），物理学家观测到了非厄米独特的层次划分：
    1. 拥有最大虚部能量（$\text{Im}(E) = 4.69$）的模式高度集中在最外层的角格点上。
    2. 随着虚部能量的降低（$\text{Im}(E) = 4.48, 4.25, 4.01$），角态波函数的最大值开始向内层晶胞的次级角格点转移（即向内侧收缩）。

结论： 张量网络算法计算出的谱线与精确对角化得到的本征值分布完全吻合，成功捕获了由于非厄米损耗诱导出的角态层次分裂物理。

2.3 十亿级超大晶格计算性能与数据（格点数 $N = 1024^3 \approx 1.07 \times 10^9$）

在证实了算法的精确性之后，研究人员将晶格线性尺寸提升至 $L_x = L_y = L_z = 1024$，体系的总自由度达到了惊人的 1,073,741,824。这是凝聚态物理计算史上首次在如此庞大的三维非厄米晶格中计算实空间局部态密度。

2.3.1 计算参数配置

  • 最大虚拟键维度： $\chi = 100$
  • 切比雪夫多项式展开阶数： $N_c = 500$
  • 空间调制： 与小体系相同的 BBH 跃迁和损耗势能函数。

2.3.2 弱损耗计算数据（$\lambda = 1.4$）

  • 复平面角态谱（图 3(a),(b)）： 在宏观极限下，通过对其中一个角晶胞（Corner Cube）的 8 个顶点进行 LDOS 求和，得到了实空间角态密度分布 $\rho_{\text{corner}}$。在复能量平面上，可清晰观测到一个明亮的孤立峰，位于 $\text{Im}(E) \approx 0.25, \text{Re}(E) = 0$ 处，表明宏观拓扑角态在体能隙中依然极其稳定。
  • LDOS 空间局域化（图 3(c)）： 展示了局部 3D 角格点胞内的电荷分布，电荷依旧高度偏置在最外角点（坐标 $(0, 0, 0)$）。

2.3.3 强损耗计算数据（$\lambda = 5$）

  • 谱线分裂现象（图 4(a),(b)）： 即使在 $1024^3$ 这样极易受到数值不稳定性和有限精度干扰的超大系统中，算法依然在复平面上极其干净地解析出了四个分裂的峰，分别位于 $\text{Im}(E) \approx 4.7, 4.5, 4.25, 4.0$ 附近，与 $16^3$ 小体系的 ED 预测结果高度一致。
  • 空间局域化模式对比（图 4(c)-(f)）：
    • 对于 $E = 4.7i$ 的模式，LDOS 完美局域在角点 $(0,0,0)$。
    • 对于 $E = 4.5i$ 的模式，波函数中心转移到相邻的 $(1,0,0), (0,1,0), (0,0,1)$。
    • 对于 $E = 4.25i$ 的模式，波函数进一步向内层格点 $(1,1,0)$ 等转移。
    • 对于 $E = 4.0i$ 的模式，电荷大部分收缩至角晶胞的内部最深格点 $(1,1,1)$。

2.3.4 计算效率与收敛性分析

由于采用了 QTCI 进行算符压缩，存储一个十亿格点的 3D 紧束缚哈密顿量 MPO 仅需要 几兆字节（MB） 的内存空间。在切比雪夫迭代过程中，通过变分截断将算符的键维度限制在 $\chi=100$，这使得计算耗时随着系统线性大小呈严格的对数级线性增长（$\mathcal{O}(\log L)$），展现了该算法在超级计算机及常规工作站上的极高实用价值。

3. 代码实现细节、复现指南与开源资源

该论文对应的开源代码已完全托管于 GitHub，使用的是 Julia 语言，并深度依赖了 Julia 生态中最强大的张量网络库 ITensors.jl。

3.1 核心算法实现步骤

复现该计算的核心在于构造 Hermitized 的 MPO，并编写支持导数算符更新的切比雪夫递归循环。以下是简化的 Julia 代码框架，展示了核心的算符生成和递推逻辑：

using ITensors

# 1. 定义二进制编码的物理自旋链基底
# 假设 L 是二进制链的长，例如 3D 1024^3 对应 3 * 10 = 30 个自旋格点
function build_binary_sites(L::Int)
  return siteinds("S=1/2", L)
end

# 2. 构造 Hermitized 哈密顿量 MPO (H~)
# 我们需要额外增加一个辅助格点放在链的头部，物理维度为 2
function build_hermitized_mpo(H_mpo::MPO, omega::ComplexF64, sites)
  L = length(sites)
  aux_site = siteind("S=1/2") # 辅助格点
  all_sites = [aux_site; sites]
  
  # 利用 ITensor 的 AutoMPO 或直接张量拼接构造 H~:
  # H~ = σ⁺ ⊗ (ωI - H) + σ⁻ ⊗ (ω*I - H⁺)
  # 在代码实现中，这等效于根据辅助格点的量子数将 H 和 H⁺ 放置在分块矩阵的非对角位置
  
  # ... (详细的张量块拼接过程) ...
  
  return H_tilde_mpo
end

# 3. 切比雪夫递归计算谱函数 (NHKPM)
# 计算 ⟨L| ∂_ω* T_{2n-1}(H~) |R⟩
function run_nhkpm_recursion(H_tilde::MPO, dH_tilde_domega::MPO, R_mps::MPS, L_mps::MPS, Nc::Int)
  # 初始化切比雪夫向量
  v_n_minus_2 = R_mps
  v_n_minus_1 = apply(H_tilde, R_mps; cutoff=1e-10, maxdim=100)
  
  # 初始化导数向量
  dv_n_minus_2 = MPS(siteinds(R_mps)) # 设为零态
  # dv_1 = ∂_ω* (H~) |R⟩
  dv_n_minus_1 = apply(dH_tilde_domega, R_mps; cutoff=1e-10, maxdim=100)
  
  spectral_sum = 0.0
  
  for n in 2:Nc
    # 递推式: v_n = 2 * H~ * v_{n-1} - v_{n-2}
    temp_v = apply(H_tilde, v_n_minus_1; cutoff=1e-10, maxdim=100)
    v_n = 2 * temp_v - v_n_minus_2
    
    # 递推式: dv_n = 2 * (∂_ω* H~) * v_{n-1} + 2 * H~ * dv_{n-1} - dv_{n-2}
    term1 = apply(dH_tilde_domega, v_n_minus_1; cutoff=1e-10, maxdim=100)
    term2 = apply(H_tilde, dv_n_minus_1; cutoff=1e-10, maxdim=100)
    dv_n = 2 * term1 + 2 * term2 - dv_n_minus_2
    
    # 如果是奇数项 (2n-1)，累加对谱函数的贡献，并乘上 Jackson 核函数因子 g_n
    if isodd(n)
      g_n = compute_jackson_kernel(n, Nc)
      term_value = inner(L_mps, dv_n)
      spectral_sum += (-1)^((n+1)÷2) * g_n * term_value
    end
    
    # 滚动更新指针
    v_n_minus_2 = v_n_minus_1
    v_n_minus_1 = v_n
    
    dv_n_minus_2 = dv_n_minus_1
    dv_n_minus_1 = dv_n
  end
  
  return (2 / pi^2) * spectral_sum
end

3.2 复现指南与开源链接

  • 官方代码仓库： GitHub - YITAOSUN42/Tensor_network_NH_tb
  • 运行环境需求：
    • Julia 1.9 及以上版本
    • 安装依赖包：Pkg.add(["ITensors", "HDF5", "LinearAlgebra", "SparseArrays"])
  • 快速复现步骤：
    1. 克隆仓库：git clone https://github.com/YITAOSUN42/Tensor_network_NH_tb.git
    2. 进入 example 目录，首先运行小尺寸基准测试：julia benchmark_small_system.jl。此脚本将计算 $16^3$ 体系的谱函数，并与内置的 ED 求解器结果进行绘图对比。
    3. 运行超大体系计算：修改配置文件中的 $L = 30$（对应三维十亿格点），执行 julia run_billion_scale.jl。该程序会自动开启多线程加速，并利用变分算符方法在数小时内输出特定实空间角点的极高精度局部谱图。

4. 关键引用文献与学术局限性批判

4.1 关键引用文献

本工作建立在非厄米物理学与张量网络算法近期多项突破的基础之上，以下五篇文献对于理解该研究的学术脉络至关重要：

  1. [1] Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitian physics, Adv. Phys. 69, 249 (2020).
    • 贡献： 提供了非厄米物理领域的权威综述，系统阐述了非厄米拓扑、异常点、增益/损耗开放系统的热力学及动力学基础。
  2. [27] G. Chen, F. Song, and J. L. Lado, Topological spin excitations in non-hermitian spin chains with a generalized kernel polynomial algorithm, Phys. Rev. Lett. 130, 100401 (2023).
    • 贡献： 首次提出了结合 Hermitization 技巧与核多项式方法（NHKPM）来求解非厄米哈密顿量谱函数的理论框架，是本工作算法的直接理论源头。
  3. [35] M. K. Ritter, Y. Núñez Fernández, M. Wallerberger, J. von Delft, H. Shinaoka, and X. Waintal, Quantics tensor cross interpolation for high-resolution parsimonious representations of multivariate functions, Phys. Rev. Lett. 132, 056501 (2024).
    • 贡献： 系统阐明了 QTCI 算法在压缩多维连续函数和实空间复杂紧束缚模型上的数学机理，使得超大晶格算符的张量化成为可能。
  4. [47] Y. Sun, M. Niedermeier, T. V. C. Antão, A. O. Fumega, and J. L. Lado, Self-consistent tensor network method for correlated super-moiré matter beyond one billion sites, Phys. Rev. Res. 7, 043288 (2025).
    • 贡献： 在厄米体系中，率先实现了超过 10 亿格点的自洽张量网络模拟，为本工作向非厄米领域的推广奠定了工程和技术基础。
  5. [51] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).
    • 贡献： 经典的核多项式方法（KPM）综述，详述了 Jackson 核、Chebyshev 递归以及实空间局部态密度计算的标准流程。

4.2 局限性批判与未来改进方向

尽管该工作在计算尺度上取得了史无前例的突破，但从学术和计算物理的角度来看，该方法仍存在以下不容忽视的局限性：

1. 强无序/无关联随机体系中的“键维度爆炸”

该算法之所以能实现十亿格点的对数级内存压缩，完全依赖于哈密顿量 $H$ 和损耗项 $V(x,y,z)$ 具有良好的低纠缠（Low-entanglement）或低多项式 rank 结构（如平移对称性、周期性调制或干净的准晶构造）。然而，如果引入真正的**强无序（Strong Random Disorder）**以研究 Anderson 定域化：

  • 实空间势能将失去任何代数关联。
  • 此时，利用 QTCI 拟合得到的 MPO 键维度 $\chi$ 将不再收敛于较小的值（如 100），而是随着系统尺寸呈指数增长。
  • 这将导致张量收缩的复杂度急剧退化，甚至可能比常规的稀疏矩阵方法还要慢。

2. 复平面二维扫描的高昂计算成本

在厄米系统中，谱函数（LDOS）仅需要在实能量轴 $\omega \in \mathbb{R}$ 上进行一维扫描。然而，非厄米系统的特征能量是复数 $\omega = \omega_r + i\omega_i$：

  • 为了绘制完整的复平面谱函数响应图，必须对复平面进行二维网格扫描（2D Grid Scan）。
  • 对于每一个复数能量点 $\omega$，都必须重新构造一次厄米化哈密顿量 $\tilde{H}(\omega)$，并重新执行完整的 $N_c = 500$ 步切比雪夫多项式递推。这使得高分辨率的二维光谱扫描极为耗时，亟需开发能“一次计算，全平面复用”的自适应插值算法。

3. 切比雪夫阶数 $N_c$ 对谱分辨率的本质限制

非厄米系统中许多最有趣的物理（如靠近异常点 EPs 时的能谱突变）表现为极其微细的谱线结构。核多项式方法的分辨率与切比雪夫展开阶数 $N_c$ 成反比（$\Delta E \sim W / N_c$，其中 $W$ 为能带总宽度）。如果要达到极高的能量分辨率，必须将 $N_c$ 提升至 $10^4$ 甚至更高：

  • 在张量网络中，高阶数的切比雪夫迭代意味着累积截断误差的急剧放大。
  • 算符的物理状态表示会逐渐失真，限制了该方法在研究极窄能隙、微扰物理或精细多体激发谱时的表现。

5. 补充讨论与跨学科展望

5.1 与量子化学及材料科学的深度融合

该张量网络算法虽然是在凝聚态紧束缚模型的语境下提出的，但其底层的“空间-二进制”编码和算符压缩哲学，对大分子量子化学模拟和宏观材料设计具有极其深远的跨学科启示价值。

1. 共价有机框架（COF）与金属有机框架（MOF）的激发态模拟

在光电、催化和储能材料中，三维共价有机框架（COF）和金属有机框架（MOF）常具有巨大的大分子晶胞，其格点数和轨道数极其庞大。在发生光吸收或电子输运时，材料不可避免地与外界环境发生能量交换（即表现为非厄米损耗和增益）。利用本文提出的 QTCI-NHKPM 框架，量子化学家可以直接构建包含数百万个原子轨道的真实 3D MOF 晶格，并在实空间直接计算特定活性位点（Active Sites）的非厄米准粒子激发表（Excitonic Spectra），而无需诉诸不切实际的周期性边界条件简化。

2. 真实多激子产生与光子传输器件设计

在现代光子学超材料中，科研人员利用非厄米性质来设计高性能的光子晶体激光器、非对称光传输器以及单向光导纤维。这些系统通常运行在经典的电磁波波动方程下，但其数学形式与紧束缚薛定谔方程完全等价。本文开发的算法能够直接帮助光电工程师在宏观尺度上设计极其复杂、具有任意几何边界的三维光子超材料，精确模拟光子流在介质中的非厄米折射、衰减与局域化行为，避免了传统时域有限差分（FDTD）在超大结构上面临的内存瓶颈。

┌─────────────────────────────────────────────────────────┐
│            QTCI-NHKPM 跨学科应用图景                    │
├──────────────────────────┬──────────────────────────────┤
│       凝聚态物理         │          量子化学            │
│  - 3D 非厄米拓扑角态     │  - 巨型分子晶体激发态        │
│  - 准晶与超莫尔超导性   │  - 活性位点非厄米准粒子谱    │
├──────────────────────────┼──────────────────────────────┤
│       光子超材料         │          声学器件            │
│  - 光子晶体激光器模拟    │  - 3D 声学拓扑绝缘体设计     │
│  - 非对称光传输器件设计  │  - 增益/损耗声波导路         │
└──────────────────────────┴──────────────────────────────┘

5.2 总结：张量网络角色的范式转变

在过去三十年中，物理学家主要将张量网络（如密度矩阵重整化群 DMRG）视为寻找**强关联量子多体系统基态（Ground State）**的专属工具。而本文所展示的突破代表了张量网络角色的关键范式转变：

张量网络正在演变为一种极具威力的大规模实空间偏微分方程（和单粒子哈密顿量）的通用数学求解器。

通过引入二进制量化编码（Quantics Representation），张量网络能够将实空间坐标转换为虚拟纠缠链，从而将空间维度转化为张量网络中的“时间”。这种变换彻底粉碎了“维度灾难”，让物理学家能够以对数级的计算成本，在宏观尺度上自由探索复杂的拓扑物态与非厄米前沿物理。我们有理由相信，随着该算法在强无序自适应、多体非厄米扩展等方向的不断演进，这一技术必将在凝聚态物理、量子信息与计算材料学领域催生出更多振奋人心的发现。