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

执行摘要

量子计算的快速演进正处于从“量子霸权”走向“实用化学与物理模拟”的过渡期。在量子化学、凝聚态物理和高能物理学中，大量的物理体系（如分子振动光谱、光子晶体、超流体以及规范场论）本质上是在无限维希尔伯特空间中展开的玻色子（Bosonic）或连续变量（Continuous-Variable, CV）系统。传统的量子计算主要基于离散变量（Discrete-Variable, DV）的量子比特（Qubit）架构，若要使用 Qubit 模拟这些玻色子系统，必须引入复杂的截断（Truncation）与映射编码（如 Gray 编码、Jordan-Wigner 类似物），这不仅伴随巨大的计算开销，还常常破坏了系统固有的代数结构。

为了克服这一根本性痛点，连续变量（CV）量子计算应运而生。它直接以谐振腔、离子振动模（Qumode）等连续物理自由度为基本信息单元。然而，如何在 CV 硬件上高效地构建变分算法以制备复杂多体系统的基态，依然是领域内的前沿难题。传统的哈密顿变分量子特征求解器（HVA-CV-VQE）由于采用固定的、受 Trotter 步长限制的硬编码线路，往往导致量子线路深度过大，极易陷入贫瘠高原（Barren Plateaus）或局部极小值。

针对这一瓶颈，最新研究提出了连续变量自适应变分量子特征求解器（CV-ADAPT-VQE）。该方法首次将自适应算子选择机制引入 CV 体系，通过动态构建高度定制化的浅层变分线路，实现了极高的基态收敛精度。本文将从核心科学问题、理论基础、技术难点、算法细节、关键基准体系数据演化，以及代码复现指南等多个维度，对这一开创性工作进行深度学术剖析。

一、 核心科学问题、理论基础与技术细节

1.1 离散变量（DV）与连续变量（CV）的范式冲突

在传统的 DV 量子计算中，两能级系统（Qubit）是基本构建块。当我们试图用 Qubit 模拟一个含有 $N$ 个玻色模（每个模最多容纳 $d$ 个玻色子）的格点系统时，所需的量子比特数呈 $\mathcal{O}(N \log_2 d)$ 增长。更严重的是，非线性的玻色相互作用算符在 Qubit 映射后会转化为复杂的、非局域的多比特 Pauli 算符乘积。这直接导致了量子线路的门复杂度极高，难以在近期的噪声中尺度量子（NISQ）设备上实现。

连续变量（CV）量子模拟则提供了一种“原生”的解决方案。在 CV 架构中，基本信息单元是Qumode（量子模），对应于一个理想的谐振子。Qumode 的状态由正则正交算符 $\hat{q}$（位置）和 $\hat{p}$（动量）描述，它们满足正则对易关系：

$$ [\hat{q}, \hat{p}] = i \mathbb{I} $$

通过玻色产生算符 $\hat{a}^\dagger$ 和湮灭算符 $\hat{a}$，我们可以定义无尽的 Fock 空间基底 $\{|n\rangle\}_{n=0}^\infty$。其动作满足：

$$ \hat{a}|n\rangle = \sqrt{n}|n-1\rangle, \quad \hat{a}^\dagger|n\rangle = \sqrt{n+1}|n+1\rangle $$

对于实际的经典数值模拟，我们通过在局部 Fock 空间中引入截断维度 $D$ 来限制计算空间，即定义 $\mathcal{H}_D = \text{span}\{|0\rangle, \dots, |D-1\rangle\}$。CV 量子模拟的核心优势在于，玻色算符是 Qumode 上的原生操作，无需进行任何映射编码即可直接对应物理相互作用。

1.2 玻色格点模型及其物理对称性

本项研究聚焦于两种最具代表性的玻色子格点模型，它们各自体现了截然不同的代数结构与对称性约束：

1.2.1 Bose-Hubbard (BH) 模型

Bose-Hubbard 模型是描述光晶格中超冷原子行为、超流-莫特绝缘体相变的核心哈密顿量。其哈密顿量可表示为：

$$ \hat{H}_{\text{BH}} = - J \sum_{r=1}^{N_S} (\hat{a}_r^\dagger \hat{a}_{r+1} + \text{h.c.}) - \frac{U}{2} \sum_{r=1}^{N_S} \hat{N}_r(\hat{N}_r - 1) $$

其中 $J$ 代表相邻格点间的跃迁（Hopping）强度，$U$ 代表格点内的玻色子排斥（$U < 0$）或吸引（$U > 0$）相互作用，$\hat{N}_r = \hat{a}_r^\dagger \hat{a}_r$ 为粒子数算符。BH 模型具有一个至关重要的物理对称性：全局 $U(1)$ 对称性，即总玻色子数 $\hat{N}_{\text{tot}} = \sum_r \hat{N}_r$ 与哈密顿量对易 $[\hat{N}_{\text{tot}}, \hat{H}_{\text{BH}}] = 0$。这意味着体系的总粒子数是严格守恒的，希尔伯特空间可以按粒子数扇区（Sector）进行精确的块对角化。

1.2.2 玻色 Kitaev 链 (Bosonic Kitaev Chain, BKC)

作为费米子 Kitaev 链的玻色子对应物，BKC 模型引入了超导配对项（或称双模挤压项），导致粒子数守恒破缺。其哈密顿量为：

$$ \hat{H}_{\text{BKC}} = - J \sum_{r} (\hat{a}_r^\dagger \hat{a}_{r+1} + \text{h.c.}) - \Delta \sum_{r} (\hat{a}_r \hat{a}_{r+1} + \text{h.c.}) - \mu \sum_{r} \hat{N}_r $$

其中 $\Delta$ 是配对相互作用强度，$\mu$ 是化学势。由于配对项 $\hat{a}_r \hat{a}_{r+1}$ 的存在，该哈密顿量不守恒总粒子数，但它保留了全局 $\mathbb{Z}_2$ 宇称对称性。宇称算符定义为 $\hat{P} = (-1)^{\hat{N}_{\text{tot}}}$，且满足 $[\hat{P}, \hat{H}_{\text{BKC}}] = 0$。基态必定落在偶宇称（Even-parity）扇区内。此外，BKC 模型是一个纯高斯哈密顿量（仅包含二次项），其基态可通过 Bogoliubov-de Gennes (BdG) 变换精确求解。

然而，当向 BKC 中引入格点内的非高斯 Kerr 相互作用时，哈密顿量演变为：

$$ \hat{H}_{\text{BKC+Kerr}} = \hat{H}_{\text{BKC}} - \frac{U}{2} \sum_{r} \hat{N}_r (\hat{N}_r - 1) $$

这一项破坏了高斯模拟的经典可积性，使得体系具有非高斯量子纠缠特性，成为检验变分算法处理复杂关联体系能力的理想试金石。

1.3 连续变量 ADAPT-VQE 算法架构

CV-ADAPT-VQE 的核心思想是摒弃固定的 ansatz 结构，通过动态、贪婪地从一个精心构造的“算子库（Operator Pool, $\mathcal{P}$）”中选择对降低能量贡献最大的算子来逐步生长（Grow）量子线路。其标准流程如下（参见论文 Fig. 1）：

  1. 初始化：制备一个参考态 $|\psi^{(0)}\rangle$。对于守恒总粒子数的 BH 模型，选择粒子均匀或交错分布的 Fock 乘积态；对于不守恒粒子数但守恒宇称的 BKC 模型，选择真空态 $|0\rangle^{\otimes N_S}$。为了避免初始梯度消失，引入“温和启动（Warm-start）”线路。

  2. 能量梯度测量：在当前迭代步 $n$，体系状态为 $|\psi^{(n)}\rangle$。针对算子库 $\mathcal{P} = \{\hat{A}_r\}$ 中的每一个生成算符 $\hat{A}_r$，评估其在零参数处的能量梯度。对于单参数算子，梯度计算公式为：

    $$ g_r = \frac{\partial E^{(n)}}{\partial \theta_r} \Big|_{\theta_r=0} = \langle\psi^{(n)}|[\hat{H}, \hat{A}_r]|\psi^{(n)}\rangle $$

    对于双参数算子（如束束器 $BS(\theta, \phi)$ 或双模挤压算子 $S_2(\theta, \phi)$），梯度采用欧几里得范数度量：

    $$ \|g_r\|_2 = \sqrt{ \left( \frac{\partial E^{(n)}}{\partial \theta_r} \right)^2 + \left( \frac{\partial E^{(n)}}{\partial \phi_r} \right)^2 } \Bigg|_{\theta_r=\phi_r=0} $$

  3. 收敛判定：若最大梯度 $\max_r \|g_r\| < \epsilon$（或能量误差满足终止条件），则算法停止。

  4. 线路生长与重参数化：选择梯度绝对值最大的算子 $\hat{O}_{n+1}$，将其插入到当前变分线路的最前端（或末端），形成新的状态表示：

    $$ |\psi^{(n+1)}\rangle = \hat{O}_{n+1}(\vartheta_{n+1})|\psi^{(n)}\rangle $$

  5. 全局参数优化：使用经典优化器（如 BFGS）对当前线路中所有的变分参数 $\boldsymbol{\vartheta}^{(n+1)} = \{\vartheta_1, \dots, \vartheta_{n+1}\}$ 进行联合优化，最小化期望值 $E^{(n+1)} = \langle\psi^{(n+1)}|\hat{H}|\psi^{(n+1)}\rangle$。之后令 $n \leftarrow n+1$，返回步骤 2。

1.3.1 对称性保护算子库（Symmetry-Preserving Pools）的设计

算子库的选择直接决定了变分搜索的效率和可行性。作者针对不同的物理对称性，专门设计了以下两套算子库（详见论文 Table I）：

  • $\mathcal{P}_{\text{BH}}$（总粒子数守恒）：仅包含粒子数守恒的算子。包括束束器 $\hat{A}_{BS, r, s} = \hat{a}_r \hat{a}_s^\dagger - \hat{a}_r^\dagger \hat{a}_s$、单模相位旋转 $\hat{A}_{R, r} = i\hat{N}_r$、单模 Kerr 门 $\hat{A}_{K, r} = i\hat{N}_r^2$ 以及双模交叉 Kerr 门 $\hat{A}_{CK, r, s} = i\hat{N}_r\hat{N}_s$。通过限制在此算子库内，量子线路绝不会跳出指定的粒子数子空间。
  • $\mathcal{P}_{\text{BKC}}$（全局宇称守恒）：包含宇称保护的 Gaussian 与非 Gaussian 算子。包括束束器、单模相位旋转、单模挤压 $\hat{A}_{S, r} = \hat{a}_r^2 - (\hat{a}_r^\dagger)^2$、双模挤压 $\hat{A}_{S_2, r, s} = \hat{a}_r\hat{a}_s - \hat{a}_r^\dagger\hat{a}_s^\dagger$。对于带有 Kerr 的模型，加入 Kerr 与交叉 Kerr 算子。

二、 关键基准体系、数据演化与性能评估

为了全面验证 CV-ADAPT-VQE 的优势，本研究通过构建基于 GPU 的高精度微分模拟器对 BH 和 BKC 模型进行了详尽的基准测试。

2.1 Bose-Hubbard 模型的数值结果

实验选择了一个 $N_S = 4$ 格点、总粒子数 $N_{ ext{tot}} = 12$ 的一维环状链（使用周期性边界条件 PBC），局部截断维度设为 $D = 13$（此时截断没有引入任何截断误差，因为单个格点上粒子数绝不可能超过 12）。

测试包含了两种极端物理物理图像（论文 Fig. 4）：

  1. 吸引相互作用极限（$\Lambda_{\text{BH}} = +4.0$）：此时粒子倾向于聚集在同一格点，基态希尔伯特空间分布较广，属于多形态纠缠结构。
  2. 排斥相互作用极限（$\Lambda_{\text{BH}} = -4.0$）：粒子倾向于均匀分布，呈现类似 Mott 绝缘体的状态，基态主要由少数几个 Fock 构型主导。
算法类型相互作用特征 $\Lambda_{\text{BH}}$达到收敛所需的门数量能量误差 $\Delta E$态不保真度 $\mathcal{I}$
HVA-CV-VQE+4.0 (吸引)160$\sim 2 \times 10^{-3}$$\sim 1 \times 10^{-3}$
CV-ADAPT-VQE+4.0 (吸引)140$\mathbf{\sim 1 \times 10^{-4}}$$\mathbf{\sim 1 \times 10^{-5}}$
HVA-CV-VQE-4.0 (排斥)160$\sim 1 \times 10^{-4}$$\sim 1 \times 10^{-4}$
CV-ADAPT-VQE-4.0 (排斥)100$\mathbf{\sim 1 \times 10^{-5}}$$\mathbf{\sim 1 \times 10^{-6}}$

核心结论：无论是在吸引还是排斥状态下，CV-ADAPT-VQE 都能以显著更少的参数和更浅的线路深度，实现比传统 HVA 高出 1 到 2 个数量级的能量精度与保真度。在 $\Lambda = -4.0$ 情况下，自适应方法在仅使用 100 个门时就达到了 HVA 使用 160 个门都无法企及的精度极限。

2.2 玻色 Kitaev 链（BKC）模型的精度优势

BKC 体系由于不守恒粒子数，局域希尔伯特空间的截断维度 $D = 6$（对应最高激发态为 $|5\rangle$）。物理参数选取为：配对强度 $\Delta = 0.5$。研究对比了两个极具代表性的区域（参见论文 Fig. 5）：

  • 不稳定性区域边界内（$\mu = -3.1$）：系统极其接近动力学不稳定边界 $\mu_c = -2(J+\Delta) = -3.0$。在此区域，基态具有强挤压（Squeezing）特性，平均粒子数 $\langle\hat{N}_{\text{tot}}\rangle$ 极大，系统关联极强。
  • 稳定区域（$\mu = -2.7$）：粒子数物理涨落较小，经典模拟极易收敛。

针对格点数 $N_S = 3, 4, 5$ 的深入评估表明：

  • 在 $N_S = 3$ 体系中：
    • CV-ADAPT-VQE 在仅仅引入约 140 个变分门 后，能量误差 $\Delta E$ 便陡降至 $10^{-10}$ 以下（逼近双精度机器极限），成功实现完美基态重构。
    • HVA-CV-VQE 则在线路增长至 150 个门时遇到了明显的收敛瓶颈，能量误差依然停留在 $10^{-4}$ 到 $10^{-3}$ 附近，发生严重平台滞留。
  • 随格点数 $N_S$ 扩展时的表现：对于 $N_S = 4$ 和 $N_S = 5$ 体系，随着量子态空间维度的指数型跃升（$5^{N_S}$），HVA 的基态重构能力迅速衰退，而不保真度在 150 个门内无法低于 $10^{-2}$。与之形成鲜明对比的是，CV-ADAPT-VQE 保持了极佳的下行通道，能量误差呈现清晰的指数衰减曲线。

2.3 BKC 与 Kerr 相互作用结合的挑战

引入非高斯 Kerr 相互作用 ($U = -1$) 后，由于纠缠态不再能被简单的辛矩阵描述，经典 Gaussian 状态变分方法彻底失效。从论文 Fig. 6 的数据来看：

在 $N_S = 3, 4, 5$ 的所有尺度下，CV-ADAPT-VQE 依旧表现出无可匹敌的鲁棒性。虽然非高斯项的引入略微拉低了前期的收敛斜率（表明非高斯基态的路径选择复杂度有所提升），但自适应算法在门数量约为 140 步时，对 $N_S=3$ 体系依然实现了 $\Delta E < 10^{-6}$ 的优异收敛。这直接确立了该算法在强关联非高斯玻色多体系统中的核心应用价值。

三、 代码实现细节、算子平铺与黑盒复现指南

本项工作的数值模拟全部是在基于 Python 的可微框架中实现的，利用了现代 GPU 加速和自动微分技术。下面对代码实现中的核心难点和优化技术进行解密。

3.1 核心代码逻辑与算子递归计算

在传统的变分模拟中，诸如挤压算子 $S(\theta, \phi) = \exp\{ \frac{\theta}{2}(e^{-i\phi}\hat{a}^2 - \text{h.c.})\}$ 的矩阵表示，通常需要对高维稀疏矩阵进行显式的矩阵指数化（Matrix Exponentiation, scipy.linalg.expm），这在算子库庞大且需要频繁求导时会造成灾难性的经典计算延迟。

为了解决这一痛点，本项工作采用了一种**稳定递推关系（Recurrence Relations）**算法。该算法不通过直接矩阵指数化，而是利用谐振子波函数在 Fock 表象下的正交多项式性质，直接递归计算出变分门在 Fock 基底下的矩阵元。其时间复杂度从 $\mathcal{O}(D^3)$ 降至 $\mathcal{O}(D^2)$，且能保持极高精度的数值稳定性。所有的算子和哈密顿量都使用稠密矩阵（Dense Matrix）在 JAX 中进行加速表达。

3.2 算子平铺（Operator Tiling）技术

当系统尺寸 $N_S$ 增大时，原始算子库的规模 $\|\mathcal{P}\|$ 也会迅速增长，这导致每一步自适应迭代中需要测量的梯度数量大幅上升。为了控制这一测量开销，作者引入了**算子平铺（Operator Tiling）**技术（见附录 B）：

  1. 小系统预演：首先在尺寸较小的经典可积系统 $N_{\text{small}}$ 上运行完整的 CV-ADAPT-VQE，记录下被高频选中的优势算子集合 $\mathcal{P}_{\text{selected}}$。
  2. 平移对称平铺：将这些优势算子沿着大系统的格点进行空间平移平铺，并用恒等算子 $\mathbb{I}$ 补齐其余格点。这相当于构建了一个专门针对该类格点模型物理特性的“自适应小算子库” $\mathcal{P}_{\text{tiled}}$。

论文 Fig. 7 证明，在 $N_S=5$ 的 BKC 系统中，相比于包含 30 个算子的完整算子库，使用平铺技术精简至 25 个算子的平铺算子库，其收敛轨迹和能量精度与完整算子库完全重合，成功实现了计算复杂度的常数级压缩。

3.3 海森矩阵回收（Hessian Recycling）

在经典的 BFGS 二阶准牛顿优化算法中，每次线路加入新算子后，优化器都需要重新估算或冷启动海森矩阵（Hessian Matrix），这消耗了大量的梯度评估步骤。作者在附录 D 中引入了海森矩阵回收（Hessian Recycling）技术：

$$ \boldsymbol{H}^{(n+1)} \approx \begin{bmatrix} \boldsymbol{H}^{(n)} & 0 \\ 0 & h_{\text{new}} \end{bmatrix} $$

将上一步优化完成的海森矩阵直接继承给下一步，作为其拟牛顿更新的初值。如论文 Fig. 10 所示，采用海森矩阵回收（HR）后，BFGS 在每一步迭代中的平均迭代次数下降了 50% 以上。这对于未来在实际量子硬件上运行该变分算法，减少物理梯度的测量开销具有决定性意义。

3.4 快速复现指南

本研究的代码已在 GitHub 开源。以下是复现核心实验的简要步骤：

环境配置与运行

# 克隆仓库
git clone https://github.com/gloriatejegar/CV_ADAPT_VQE.git
cd CV_ADAPT_VQE

# 创建虚拟环境并安装依赖
conda create -n cv_adapt python=3.10 -y
conda activate cv_adapt
pip install --upgrade "jax[cuda11_pip]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html
pip install scipy numpy

# 运行 Bose-Hubbard 模型自适应优化模拟 (示例命令行参数)
python run_adapt_vqe.py --model bose_hubbard --nsites 4 --nparticles 12 --cutoff 13 --lambda_bh 4.0

代码结构解析

  • simulator/：包含基于 JAX 的原生 CV 门控模拟引擎，内部实现了递归计算矩阵元的核心逻辑。
  • pools.py：定义了各类对称性保护算子库，包括 $\mathcal{P}_{\text{BH}}$ 和 $\mathcal{P}_{\text{BKC}}$ 的算符生成器。
  • adapt.py：包含自适应选择主循环、梯度评估（Euclidean norm 梯度计算）以及终止条件判定。
  • optimizers.py：集成了 BFGS 优化器，并支持温和启动（Warm-start）及海森矩阵回收（Hessian Recycling）机制。

四、 学术局限性批判、关键文献追溯与未来展望

4.1 本项工作的局限性分析

尽管 CV-ADAPT-VQE 取得了令人瞩目的基准测试结果，但作为一项开创性探索，它在走向实际大规模应用前仍存在以下不容忽视的局限性：

  1. 局部 Fock 截断误差与幺正性丧失：在对不守恒粒子数的哈密顿量（如 BKC 系列）进行模拟时，人为截断 $D$ 必然会导致一些高阶玻色激发态被强行抹去。当应用算子库中那些能够改变粒子数的算子（如挤压门）时，在截断边界处的矩阵元素会发生溢出，从而破坏算符的严格幺正性（Unitarity）。虽然作者在每步门操作后通过人工重新归一化状态（State Normalization）进行了对冲，但这只是数值上的妥协，在真实量子硬件上这种幺正性的偏离可能导致状态制备失效。
  2. 经典稠密矩阵模拟器的内存瓶颈：由于 JAX 内部缺乏对超大规模高阶张量稀疏矩阵的自动微分原生支持，代码采用了稠密矩阵表象。这导致内存开销以 $\mathcal{O}(D^{2N_S})$ 极速暴涨。目前最大模拟尺寸被锁死在 $N_S = 5, D = 6$。若要进一步研究凝聚态物理中的热力学极限（Thermodynamic Limit），必须开发基于矩阵乘积态（MPS）或稀疏张量网络的可微 CV 模拟器。
  3. 非高斯温和启动的局部陷阱依赖性：由于在 Fock 基底下的对角算子（如 Kerr 项）在任何 Fock 乘积态上的初始梯度都严格为零，CV-ADAPT-VQE 面临着极强的“冷启动失效”风险。虽然本工作通过手动添加一层 HVA 作为温和启动（Warm-start）成功规避了这一问题，但这种温和启动的设计高度依赖对目标哈密顿量基态先验物理性质的理解。如果面对的是一个完全未知的量子化学反应分子势能面，如何设计无偏置（Unbiased）且高效的启动机制依然是一个悬而未决的难题。

4.2 关键文献关联与理论源流

要深刻理解本项工作，以下四篇里程碑式的文献构成了其理论支柱：

  • 文献 [9] (S. Lloyd & S. L. Braunstein, Phys. Rev. Lett. 1999)：奠定了连续变量通用量子计算的理论基石，首次证明了高斯门配合单个非高斯算子（如立方相位门或 Kerr 门）即可实现通用的 CV 量子计算。
  • 文献 [56] (H. R. Grimsley et al., Nat. Commun. 2019)：首次提出了自适应变分量子特征求解器（ADAPT-VQE）框架，开创了动态生长低深度变分量子线路的先河。
  • 文献 [68] (S. Yalouz et al., Quantum 2021)：首次尝试将哈密顿变分拟设（HVA）引入 CV 量子模拟中并应用于 Bose-Hubbard 模型，这也是本项工作直接对标和超越的硬编码基准。
  • 文献 [101] (J. S. Van Dyke et al., Phys. Rev. Res. 2024)：提出了算子平铺（Operator Tiling）技术以降低变分量子算法在大尺度系统上的测量和寻址复杂度，为本文的附录 B 奠定了理论基础。

4.3 量子化学中的未来展望：振动 VQE (vVQE)

在量子化学中，分子的电子结构通常用费米子建模。然而，分子振动光谱（Vibrational Structure）、非谐振动势能面等问题，在本质上是多模耦合非谐振子系统的运动，是典型的玻色子体系。传统的量子化学计算（如 VSCF、VCI、VCC）在处理强非谐性分子（如水团簇、大分子生物酶活性中心）时，会面临维数灾难。

CV-ADAPT-VQE 为解决分子振动问题提供了完美的范式： 每一个分子的正则振动模（Normal Mode）可以直接映射为一个硬件原生的 Qumode。通过精心构造包含非谐性校正项（如立方、四次方非谐性算符）的非高斯算子库 $\mathcal{P}_{\text{vib}}$，CV-ADAPT-VQE 能够自动、自适应地识别分子振动中的关键非谐性耦合通道，构建出极度精简的振动变分态。这不仅能极大地提高分子振动能级和红外/拉曼光谱的模拟精度，还彻底免去了繁琐的费米子-量子比特编码。该方法有望成为未来量子化学模拟强关联振动多体问题的新一代标准工具。

五、 补充：连续变量量子门代数理论与物理边界严格证明

为了使本文对技术背景深厚的研究人员更具参考价值，下面提供连续变量变分量子门控的核心数学代数表示，以及玻色 Kitaev 链稳定性边界的严格 Bogoliubov-de Gennes 理论证明。

5.1 变分算子门代数理论表示

在 CV-ADAPT-VQE 算子库中，各变分门控的生成算符和么正变换矩阵的代数形式定义如下：

1. 相位旋转门（Phase Rotation Gate）

$$ \hat{R}(\theta) = \exp\{ i \theta \hat{N} \} $$

在 Fock 基底下为对角阵：$\hat{R}(\theta)|n\rangle = e^{i \theta n}|n\rangle$。

2. 位移门（Displacement Gate）

$$ \hat{D}(\alpha) = \exp\{ \alpha \hat{a}^\dagger - \alpha^* \hat{a} \} $$

其中 $\alpha = \theta e^{i \phi}$，其将真空态 $|0\rangle$ 转化为相干态 $|\alpha\rangle$。

3. 单模挤压门（Single-Mode Squeezing Gate）

$$ \hat{S}(z) = \exp\{ \frac{1}{2}(z^* \hat{a}^2 - z (\hat{a}^\dagger)^2) \} $$

其中 $z = \theta e^{i \phi}$。它引入了单个 Qumode 位置与动量涨落的非对称压缩。

4. 束束器门（Beam Splitter Gate）

$$ \hat{BS}(\theta, \phi) = \exp\{ \theta (e^{i \phi} \hat{a}_1 \hat{a}_2^\dagger - e^{-i \phi} \hat{a}_1^\dagger \hat{a}_2) \} $$

这是双模粒子数守恒的相干混模操作，在 Fock 空间内具有块对角结构。

5. 双模挤压门（Two-Mode Squeezing Gate）

$$ \hat{S}_2(\theta, \phi) = \exp\{ \theta (e^{i \phi} \hat{a}_1 \hat{a}_2 - e^{-i \phi} \hat{a}_1^\dagger \hat{a}_2^\dagger) \} $$

产生双模间的 EPR 纠缠，同时严格保持系统总粒子数的宇称。

6. Kerr 门与交叉 Kerr 门（Non-Gaussian Gates）

$$ \hat{K}(\theta) = \exp\{ i \theta \hat{N}^2 \}, \quad \hat{CK}(\theta) = \exp\{ i \theta \hat{N}_1 \hat{N}_2 \} $$

提供非线性相位移动，是产生非高斯量子魔力（Quantum Magic）的核心源泉。

5.2 玻色 Kitaev 链不稳定性边界的严格数学证明

玻色 Kitaev 链哈密顿量在引入配对相互作用后，若参数选择不当，其能谱可能变得无界，导致系统出现动力学不稳定性（即玻色子数无限凝聚）。以下通过 Bogoliubov-de Gennes (BdG) 形式给出严格的物理稳定性边界推导：

设一维 BKC 哈密顿量为：

$$ \hat{H}_{\text{BKC}} = - J \sum_{r} (\hat{a}_r^\dagger \hat{a}_{r+1} + \hat{a}_{r+1}^\dagger \hat{a}_r) - \Delta \sum_{r} (\hat{a}_r \hat{a}_{r+1} + \hat{a}_{r+1}^\dagger \hat{a}_r^\dagger) - \mu \sum_{r} \hat{a}_r^\dagger \hat{a}_r $$

引入空间傅里叶变换：$\hat{a}_r = \frac{1}{\sqrt{N_S}} \sum_k e^{i k r} \hat{a}_k$，可将哈密顿量写为动量空间的对角块形式：

$$ \hat{H}_{\text{BKC}} = \frac{1}{2} \sum_k \begin{pmatrix} \hat{a}_k^\dagger & \hat{a}_{-k} \end{pmatrix} \mathcal{H}_{\text{BdG}}(k) \begin{pmatrix} \hat{a}_k \\ \hat{a}_{-k}^\dagger \end{pmatrix} - \frac{1}{2} N_S \mu $$

其中，BdG 矩阵为：

$$ \mathcal{H}_{\text{BdG}}(k) = \begin{pmatrix} \xi_k & \Delta_k \\ \Delta_k & \xi_k \end{pmatrix} $$

此处定义：

$$ \xi_k = -2J \cos k - \mu, \quad \Delta_k = -2\Delta \cos k $$

为了求得准粒子激发的能谱，我们需要通过辛么正变换（Symplectic Transformation）对具有非对易关系的玻色产生/湮灭算符进行对角化。对角化条件对应于求解广义特征值问题：

$$ \det \left( \sigma_z \mathcal{H}_{\text{BdG}}(k) - E(k) \mathbb{I} \right) = 0 $$

其中 $\sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$ 为第三个 Pauli 矩阵。代入计算：

$$ \sigma_z \mathcal{H}_{\text{BdG}}(k) = \begin{pmatrix} \xi_k & \Delta_k \\ -\Delta_k & -\xi_k \end{pmatrix} $$

其特征值方程展开为：

$$ \det \begin{pmatrix} \xi_k - E(k) & \Delta_k \\ -\Delta_k & -\xi_k - E(k) \end{pmatrix} = 0 \implies -(\xi_k^2 - E(k)^2) + \Delta_k^2 = 0 $$

解得准粒子激发的能谱为：

$$ E(k) = \pm \sqrt{\xi_k^2 - \Delta_k^2} = \pm \sqrt{(-2J \cos k - \mu)^2 - (2\Delta \cos k)^2} $$

稳定性分析物理判据： 系统的激发电能 $E(k)$ 必须是严格的实数。若 $E(k)$ 出现虚数，则意味着系统的时间演化算符 $\exp\{-i \hat{H} t\}$ 中将包含指数增长项 $\exp\{|\text{Im}(E)| t\}$，导致格点上的玻色子数目发生雪崩式无限增长，即发生物理系统坍塌（不稳定相变）。

因此，系统保持动力学稳定的充要条件是：

$$ \xi_k^2 > \Delta_k^2 \implies \xi_k > |\Delta_k| $$

代入 $\xi_k$ 与 $\Delta_k$ 的定义式：

$$ -2J \cos k - \mu > |2\Delta \cos k| $$

此不等式必须对布里渊区内所有的动量 $k \in [-\pi, \pi]$ 同时恒成立。为了找出最严苛的不稳定临界点，我们考察 $\cos k = 1$（即 $k=0$ 的 $\Gamma$ 点，或 $\cos k = -1$ 的边缘点）。

在 $k=0$ 处：

$$ -2J - \mu > 2\Delta \implies \mu < -2(J + \Delta) $$

这一结论与论文中给出的物理稳定性边界公式 (17) 完全契合。这一严格推导不仅指明了变分算法在探索临界相变点附近的物理约束，也深刻解释了图 2 中所展示的物理量平均激发粒子数 $\langle\hat{N}_{\text{tot}}\rangle$ 随参数跃迁而急剧暴涨的深层物理诱因。