来源论文: https://arxiv.org/abs/2606.20313v1 生成时间: Jun 19, 2026 01:28

深度解析亚欧姆自旋-玻色模型的动力学相与纠缠结构：基于树张量网络态（TTNS）与TDVP-PS的研究

0. 执行摘要

自旋-玻色模型（Spin-Boson Model, SBM）是研究开放量子系统耗散动力学、多体物理及量子退相干的核心范式。特别是在亚欧姆（sub-Ohmic, $0 < s < 1$）机制下，低频环境模式极强的耦合强度与极慢的响应时间，引发了非平凡的非平衡相变与长记忆效应。传统上，该领域的动力学相分类（相干、非相干、伪相干）主要建立在自旋极化布居数 $\langle\hat{\sigma}_z(t)\rangle$ 的时间演化特征上，然而，这些动力学相是否对应于本质上不同的自旋-环境纠缠结构，长期以来一直缺乏微观且全面的量化回答。

本研究利用先进的树张量网络态方法结合投影分裂随时间变化变分原理（TTN-TDVP-PS），对整个亚欧姆 $(s, \alpha)$ 参数空间进行了系统性的扫描。研究发现，自旋纠缠熵 $S_{\text{spin}}(t)$ 能够以远快于极化弛豫的速度快速稳定在平稳高原值 $S_{\text{stable}}$。利用这一发现构建的“静态纠缠熵景观图”显示：在低 $s$ 区域，纠缠熵最大值脊线（Entropy Ridge）与基于布居数的动力学边界高度吻合；但在大 $s$ 区域，布居数表现出的双支分叉结构（相干-非相干、非相干-伪相干）在纠缠熵景观中并未重现，纠缠熵脊线依然保持单值，并深嵌于非相干区内部。通过 Bloch 球几何表示，三种动力学相表现出独特的轨迹形态（螺旋缠绕、单调趋近及钩状受阻激发）。此外，模式解析环境纠缠分析定量证实了低频模式的主导地位，并首次揭示了相干动力学对环境模间关联（Bath-Mode Correlations）的主动增强效应。本工作不仅拓宽了利用纠缠测量表征耗散量子相变的视角，也为大自由度复杂环境量子动力学模拟树立了新的方法论标杆。

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

1.1 核心科学问题：布居数分类与纠缠结构的微观映射

在开放量子动力学中，系统与环境的相互作用会导致局域退相干和能量耗散。对于经典的两能级系统耦合谐振子库（即自旋-玻色模型），前人（如 Otterpohl 等人，Phys. Rev. Lett. 2022）通过自旋极化差 $\langle\hat{\sigma}_z(t)\rangle$ 的演化特征，将其动力学划分为三大区间：

相干区间（Coherent Regime）：表现为衰减的简谐振荡，在零点上下多次穿梭，对应于局域态间的相干隧穿。
非相干区间（Incoherent Regime）：表现为单调衰减，无任何振荡，说明强耗散已完全抑制了相干隧穿。
伪相干区间（Pseudo-coherent Regime）：在极短时间内出现一个单一的极小值，随后系统出现部分极化恢复并迅速锁死在局域态，属于强耦合诱导的极快速过阻尼局域化。

然而，这种分类完全依赖于自旋（子系统）的单体可观测属性。量子信息科学指出，子系统的退相干本质上是其与环境谐振子发生纠缠（产生不可逆关联）的结果。这就引出了核心科学问题：自旋与环境（浴）之间的纠缠熵及纠缠结构，在微观上是否与上述三大动力学相存在一一映射关系？这种纠缠在时域上如何建立？在参数空间中其静态景观呈现何种拓扑特征？

1.2 理论基础：亚欧姆自旋-玻色模型

自旋-玻色模型的总哈密顿量可写为：

$$\hat{H}_{\text{SB}} = -\frac{1}{2}\hbar \Delta \hat{\sigma}_x + \frac{1}{2}\epsilon \hat{\sigma}_z + \sum_j \left( \frac{1}{2}m_j \omega_j^2 \hat{x}_j^2 + \frac{\hat{p}_j^2}{2m_j} \right) + \frac{1}{2}\hat{\sigma}_z \sum_j c_j \hat{x}_j$$

其中 $\Delta$ 为自旋的隧穿矩阵元，$\epsilon$ 为能级偏置（本文主要研究无偏置系统 $\epsilon=0$），$\hat{\sigma}_x, \hat{\sigma}_z$ 为 Pauli 算符。第三项为环境谐振子集合，第四项为自旋与环境的线性位移耦合，耦合强度由常数 $c_j$ 控制。

系统动力学和物理行为完全由其谱密度函数（Spectral Density Function）决定：

$$J(\omega) = \sum_j \frac{c_j^2}{2\omega_j} \delta(\omega - \omega_j)$$

在连续介质极限下，谱密度通常采用幂指数截止形式：

$$J(\omega) = 2 \alpha \frac{\omega^s}{\omega_c^{s-1}} e^{-\omega/\omega_c}$$

其中 $\alpha$ 为无量纲自旋-环境耦合强度，$s$ 为谱指数，$\omega_c$ 为截止频率。根据 $s$ 的取值，模型分为：

$s > 1$：超欧姆（Super-Ohmic）机制；
$s = 1$：欧姆（Ohmic）机制；
$0 < s < 1$：亚欧姆（Sub-Ohmic）机制。

在亚欧姆机制下，谱密度在低频极限（$\omega \to 0$）处发散，这意味着系统的物理行为极大地受控于极低频的红外发散浴模式。由于这些低频模式对应着极长的时间尺度，系统表现出极强的非马尔可夫（Non-Markovian）记忆效应，这使得传统的摄动理论、马尔可夫近似甚至是传统的数值传播方法全部失效。

1.3 技术难点与传统方法的瓶颈

数值精确求解亚欧姆自旋-玻色模型在整个 $(s, \alpha)$ 空间下的长时动力学是一项公认的技术挑战：

路径积分方法（如 QUAPI）：虽然可以严格处理耗散，但在红外发散（低 $s$）时，环境关联时间趋于无限长，导致计算所需的存储器记忆时间（Memory Time Window）急剧膨胀，计算复杂度随时间步数呈指数增长。
层级运动方程（HEOM）：需要对谱密度进行有理逼近，在低温及亚欧姆发散下需要极深、极其庞大的层级结构，导致基底数量爆炸。
多态多组态随时间变化哈特里方法（ML-MCTDH）：虽然利用树状剪枝在一定程度上缓解了维数灾难，但在处理低频模式的长时振荡时，活性基矢选择和变分运动方程的奇异性（Singularity of Equations of Motion）经常导致数值积分器步长骤减甚至崩溃。
一维张量网络（矩阵乘积态，MPS/DMRG）：将连续浴离散化后，若采用一维链状网络（Star-to-Chain transformation），自旋与浴的远程耦合转化为一维链上的近邻耦合，但此变换会产生长程纠缠积累，导致 MPS 的键维数（Bond Dimension）在长时演化中迅速饱和甚至爆炸，且无法直观解析浴内模间的纠缠分布。

1.4 方法细节：树张量网络态（TTNS）与 TDVP-PS 算法

为了克服上述瓶颈，本工作引入了树张量网络态（TTNS）搭配投影分裂随时间变化变分原理（TDVP-PS）。

1.4.1 树张量网络拓扑优势

相比于 MPS 链式结构，TTNS 允许网络在拓扑上呈现层级分支树状。对于自旋-玻色模型，自旋或最重要的慢浴模式可以置于树根（Root）或靠近中心的分支节点，而其余的 $N_b = 1000$ 个离散玻色子模式则按照频率等级或物理关联度，被逐级排布在树的叶子节点（Leaves）上（其物理排布示意图可参考层次化结构树）。这种多叉树状网络结构不仅大大缩短了任意两个谐振子模式之间的网络路径距离（Path Distance），而且能够极其自然地捕获局部浴模关联与系统-浴的多体纠缠，从而维持极低的虚拟键维数。此外，网络中引入了纯虚拟的辅助节点（Auxiliary Nodes），极大地增强了高维波函数的张量分解效率。

1.4.2 自动 TTNO 构建（Bipartite-Graph Approach）

由于树拓扑分支多样，手动构建多体算符的树张量网络算符（TTNO）极其繁琐且极易出错。本研究采用基于二分图（Bipartite-graph）的自动构建算法，将哈密顿量中所有加和形式的单体和双体算符项自动、精确且极其紧凑地转化为相应的 TTNO 表示，彻底解决了变分传播中的算符作用瓶颈。

1.4.3 投影分裂时间步传播（TDVP-PS）

量子态 $\|\Psi(t)\rangle$ 的严格时间演化受控于薛定谔方程。TDVP 将其投影至 TTNS 所定义的高维非线性变分流形 $\mathcal{M}$ 的切空间上：

$$i\hbar \frac{\partial}{\partial t} \|\Psi(t)\rangle = \hat{P}_{T_{\Psi}\mathcal{M}} \hat{H}_{\text{SB}} \|\Psi(t)\rangle$$

其中 $\hat{P}_{T_{\Psi}\mathcal{M}}$ 为流形在 $\|\Psi(t)\rangle$ 处的切空间投影算符。对于树状张量网络，其投影算符可表示为一系列作用在节点（Nodes）和连接键（Edges）上的局部投影算符的交替相减形式：

$$\hat{P}_{T_{\Psi}\mathcal{M}} = \sum_{k \in \text{Nodes}} \hat{P}_k - \sum_{e \in \text{Edges}} \hat{P}_e$$

**投影分裂方法（Projector-Splitting, PS）**将整体投影时间演化算符 $e^{-i\hat{H}\delta t}$ 分解为在各节点张量上局部演化算符的乘积顺序。在一轮积分内：

在正向传播中，对某个局部节点进行局部哈密顿量对角化传播（向前推进 $\delta t$）；
紧接着在连接键上进行向后演化（向后推进 $-\delta t$），以便为下一个节点的传播更新提供准确的物理基矢。

这种分裂方案具有多重无与伦比的数学与物理优势：

无条件数值稳定性（Unconditional Stability）：即使系统包含极宽阶跃的能量范围（如自旋能级与红外/紫外谐振子能量相差数个数量级），依然能保持极其稳健的演化；
辛对称与幺正保持（Unitary & Symplectic Preserving）：严格保持波函数的归一化，且在传播过程中没有人工虚部耗散；
动态适应键维数：通过对分裂步中奇异值分解（SVD）的动态截断，能够根据设定的截断精度（$\chi_{\text{err}}$）自适应地调整树网络中每条虚拟键的物理维度。

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

2.1 模拟参数设置与物理离散化

本研究选取的系统参数完全对标经典的亚欧姆自旋-玻色模型 Benchmark 标准：

自旋隧穿强度：$\Delta = 1$（作为所有能量和时间标度的归一化基准）；
能级偏置：$\epsilon = 0$（对称双势阱）；
环境谱密度截止频率：$\omega_c = 10 \Delta$（强宽带浴限制）；
物理谐振子个数：$N_b = 1000$。这是一个极高保真度的离散化方案，能确保在较长的时间演化内不发生人工庞加莱复归（Poincaré Recurrence）。

为了极其精确地捕捉亚欧姆低频发散，谱密度离散化采用对数-线性混合密度分布（密度函数 $\rho(\omega)$）：

$$\rho(\omega) = \frac{N_b + 1}{\omega_c} e^{-\omega/\omega_c}$$

通过对累积频数函数进行数值积分 $\int_0^{\omega_j} d\omega \rho(\omega) = j$（$j=1,\dots,N_b$）生成各离散模式的频率 $\omega_j$。离散模式与自旋的耦合强度由下式决定：

$$c_j^2 = \frac{2}{\pi} \omega_j \frac{J(\omega_j)}{\rho(\omega_j)}$$

离散化的有限截止效应通过在变分哈密顿量中引入绝热重整化隧穿元 $\Delta_{\text{eff}}$ 予以严格修正：

$$\Delta_{\text{eff}} = \Delta \exp \left[ -\int_{\omega_{\text{max}}}^{\infty} d\omega \frac{J(\omega)}{\omega^2} \right]$$

其中 $\omega_{\text{max}} = \omega_{N_b}$ 为截断的最大频数。

在 TTNS 中，每个谐振子的局部玻色子基底维度限制为 $d = 16$，在强耦合区域增加到 $d=32$，通过截断检查确保其占有数几率（Occupancy Probability）低于 $10^{-8}$，防止发生人工基底饱和效应。

2.2 极化动力学与相图的复现（对标经典）

通过高精度时间演化，首先复现了自旋极化布居数 $\langle\hat{\sigma}_z(t)\rangle = \langle\Psi(t)\| \hat{\sigma}_z \|\Psi(t)\rangle$。如图 1 所示，通过在整个 $(s, \alpha)$ 面上扫描：

在固定 $s=0.7$ 时，随着 $\alpha$ 从 $0.1$ 增大到 $0.8$，系统动力学表现出明显的从相干阻尼振荡（$\alpha=0.10$）到非相干单调衰减（$\alpha=0.40$），再到伪相干受阻快速局域化（$\alpha=0.80$）的跃迁。
固定耦合强度 $\alpha=0.80$，当谱指数 $s$ 从 $0.1$ 增加到 $0.7$ 时，原本在 $s=0.1$ 强非相干锁死的状态，在 $s=0.7$ 下逐渐恢复了一定程度的隧穿，这对应于高谱指数下红外物理权重的弱化。
绘制出的动力学相图（图 1e）与经典精确 QUAPI 方法的结果（Otterpohl 等人，2022）在边界线位置高度重合，有力验证了 TTN-TDVP-PS 算法在数值上的极高精度。研究指出，在长时传播（$t = 50$ a.u.）下，TTNS 的高稳健性甚至修正了部分极窄相边界的微小数值偏离。

2.3 自旋纠缠熵动力学及“静态熵景观图”的非平凡发现

利用 TTNS 波函数，通过部分跟踪浴自由度，直接计算自旋子系统的减缩密度矩阵：

$$\rho_{\text{spin}}(t) = \text{Tr}_{\text{bath}} \|\Psi(t)\rangle\langle\Psi(t)\|$$

随时间演化的自旋纠缠熵（冯·诺依曼熵）定义为：

$$S_{\text{spin}}(t) = -\text{Tr} \left[ \rho_{\text{spin}}(t) \log \rho_{\text{spin}}(t) \right]$$

研究结果呈现出以下核心物理发现：

2.3.1 极快的高原建立速度（尺度分离现象）

在几乎整个参数范围内，纠缠熵 $S_{\text{spin}}(t)$ 在极短的时间内（$t \approx 1 \sim 3$ a.u.）便迅速攀升并锁定在一个长时间极其稳定的常数高原（Stationary Plateau）上，而此时自旋的极化动力学 $\langle\hat{\sigma}_z(t)\rangle$ 仍在大范围地剧烈振荡（如 $s=0.7, \alpha=0.1$ 的相干振荡可一直持续至 $t > 20$ a.u.）。这一极强的时间尺度分离现象表明：自旋与外部环境的多体关联早在局域布居数实现动力学平稳前，就已经达到了饱和状态。自旋之后的长时隧穿运动实际上发生在一张已经“稳固建立”的固定系统-浴纠缠网中。

2.3.2 静态纠缠熵景观图（Stationary Entropy Landscape）与分叉缺失

通过对平稳期（$t = 10 \sim 50$ a.u.）的纠缠熵进行时间平均：

$$S_{\text{stable}}(s, \alpha) = \frac{1}{N_{\text{avg}}} \sum_{t_n \ge t_{\text{avg}}} S_{\text{spin}}(t_n)$$

研究者在 $(s, \alpha)$ 参数面上绘制出了极其精美的二维纠缠熵标量场（图 2e）。该图展现了前所未有的物理特征：

熵脊线与边界的耦合：在小 $s$ 区域（$s < 0.4$），静态纠缠熵的最大值脊线 $s_{\text{ridge}}(\alpha) = \text{arg\,max}_s S_{\text{stable}}(s,\alpha)$ 几乎与布居数定义的相干-非相干相边界完全吻合。
大 $s$ 区域的分叉缺失（核心发现）：在 $s > 0.4$ 的区域，布居数动力学分裂为两支分叉边界（“相干-非相干”与“非相干-伪相干”），然而静态纠缠熵景观的脊线并未分叉。这条唯一的纠缠熵脊线依然呈单值形态，直接从非相干区域内部穿过，而没有去追踪这两支边界。这说明大 $s$ 下，虽然布居数的演化形态呈现多重分叉，但自旋-浴的底层多体关联并未在边界处发生本质的结构性突变，纠缠最强的极值依然由于低频模式权重的淡化而收敛在非相干区中央。
马鞍区 dip（图 2f）：在 $s \approx 0.4$ 附近存在一个明显的纠缠熵局部凹陷（dip）。这是由于该区域自旋与多频浴模式的耦合竞争，产生了极长的弛豫与热化时间尺度，导致在当前物理观测时间窗（$t=50$ a.u.）内，多体状态还处于一种复杂的亚稳态中。

2.4 Bloch 球动力学的几何学诠释

量子两能级系统的任意状态均可在 Bloch 球（Bloch Sphere）空间中几何表征。其减缩密度矩阵表述为：

$$\rho_{\text{spin}}(t) = \frac{1}{2} \left( \mathbf{1} + \langle\hat{\sigma}_x(t)\rangle \hat{\sigma}_x + \langle\hat{\sigma}_y(t)\rangle \hat{\sigma}_y + \langle\hat{\sigma}_z(t)\rangle \hat{\sigma}_z \right)$$

其对应的 Bloch 向量模长，即 Bloch 半径为：

$$r(t) = \sqrt{\langle\hat{\sigma}_x(t)\rangle^2 + \langle\hat{\sigma}_y(t)\rangle^2 + \langle\hat{\sigma}_z(t)\rangle^2}$$

自旋纠缠熵是 $r$ 的单调递减函数：当自旋处于纯态时，$r=1, S_{\text{spin}}=0$；当自旋与环境完全纠缠达到最大混态时，$r=0, S_{\text{spin}} = \ln 2 \approx 0.693$。静态纠缠熵进入平稳高原这一物理发现，在几何上对应于：Bloch 轨迹最终完全受限在某一个固定半径的同心球壳上运动。这一过程通过无偏置条件下末态隧穿电流 $\langle\hat{\sigma}_y(t)\rangle \to 0$ 得到了定量佐证，说明末态轨迹严格收敛在 $xz$ 平面的圆弧上。

几何轨迹分析（图 3 与图 4）对三种动力学相给出了极具穿透力的物理解释：

相干轨迹（Coherent）：Bloch 向量在三维空间中呈现出明显的绕心旋转、螺线缠绕（Spiral）形态，在经历多次螺旋收缩后，最终温和地收敛并静止在较外层的球壳上（$r \approx 0.50 \sim 0.63$，中等强度纠缠）。
峰值纠缠非相干轨迹（Peak Incoherent）：轨迹以最短路径直奔 Bloch 球心而去，最终定位在半径最小的球壳上（最小半径 $r_{\text{min}} \approx 0.39 \sim 0.48$），标志着自旋与环境建立了最强大的量子纠缠关联。
伪相干轨迹（Pseudo-coherent）：在以往的极化观察中，该相仅表现为一个极简单的单调下降。而三维 Bloch 几何（图 4d）首次揭示了其隐藏的“秘密”：Bloch 向量在初始时展现出一个向负 $x$ 方向巨大甩出的隧穿意图（$\langle\hat{\sigma}_x\rangle \approx -0.59$），但其瞬间被强悍的浴模式所锚定并生硬地“拉回”，在 $xz$ 投影面上画出一个极具代表性的钩状轨迹（Hook-shaped failed tunneling excursion）。这在几何上直观印证了伪相干本质上是一种“受阻的隧道激发”。

2.5 模式解析纠缠分析：模间关联的非平凡中介效应

为进一步透视“浴”内部的微观演化，本工作通过部分跟踪其余模式，计算了单个谐振子模 $j$ 的单模纠缠熵 $S_j(t)$，以及自旋与该浴模的双体纠缠熵 $S_{\text{spin},j}(t)$，进而构建了自旋-浴模互信息（Spin-Boson-Mode Mutual Information）：

$$I_{\text{spin}:j}(t) = S_{\text{spin}}(t) + S_j(t) - S_{\text{spin},j}(t)$$

2.5.1 $S > I$ 的物理启示与模间关联

对于单个谐振子模式，$S_j(t)$ 代表它与系统及其他所有浴模组成的“其余部分”的整体量子纠缠；而 $I_{\text{spin}:j}(t)$ 则仅仅衡量它与自旋之间的直接关联。通过分析图 5 发现：

在非相干与伪相干相中：$S_j(t) \approx I_{\text{spin}:j}(t)$。这说明在强耗散条件下，浴模的纠缠几乎完全来自于它同自旋的直接相互作用，浴模之间不存在间接量子纠缠。
在相干相中（核心物理发现）：$S_j(t) > I_{\text{spin}:j}(t)$。这表明单模总纠缠显著超越了其与自旋直接纠缠的上限。由于哈密顿量中完全不包含任何浴模之间的直接耦合项（Bath-bath coupling），该差值 $\Delta S = S_j - I_{\text{spin}:j}$ 只能是由相干演化的自旋作为物理中介，间接在不同环境模式之间介导并建立的“模间多体量子关联”（mediated bath-mode correlations）。这是首次在亚欧姆自旋-玻色模型中定量观测到相干自旋运动对环境内部纠缠网络的反作用构建功能。

2.5.2 频域标度律（Frequency Scaling）与低频主导定量化

如图 6 所示，谐振子模纠缠振荡频率 $\Omega$ 严格满足 bare 标度率：$\Omega \approx \omega_{\text{boson}}/(2\pi)$，证明其振荡完全受其裸谐振频率控制，不受极化重整化影响。而其振荡幅度 $S_{\text{Amp}}$ 随着频率 $\omega_{\text{boson}}$ 的升高而呈现极其剧烈的指数级衰减。这从微观时域演化上无懈可击地定量证实了：亚欧姆自旋-环境纠缠的绝对尺度，完全被极少数低频慢浴模式所决定，它们是退相干和耗散的根本源头。

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

3.1 所用的关键开源软件包

本研究所使用的全部数值模拟和算法流程，均基于开源张量网络 Python 框架 Renormalizer。该软件包专为分子光谱、多体纠缠和开放量子系统动力学设计，高度集成了 TTNS 拓扑、自动 TTNO 算符图论构建以及 GPU 硬件加速算法。

项目 GitHub 仓库: Renormalizer Development Group
核心功能: 基于自动多重二分图构建 MPO/TTNO；高稳定性、支持矩阵乘积态与树张量网络的 1-site & 2-site TDVP 算法。

3.2 模拟复现 Python 脚本实例

以下提供一段高度结构化的 Python 脚本，演示如何利用 Renormalizer 从头构建亚欧姆谱密度、离散化 $N_b=1000$ 个模式、构建 TTNS 树拓扑并运行高精度的 TDVP 动力学传播，最终提取自旋纠缠熵。

import numpy as np
from renormalizer.model import Model, Mol, Phonon
from renormalizer.model.basis import BasisHalfSpin, BasisOscillator
from renormalizer.mps import Mps, Mpo, Ttns, Ttno, TdvpPS
from renormalizer.utils import Op, Constant

def get_sub_ohmic_spectral_density(w, s, alpha, w_c):
    # J(w) = 2 * alpha * w^s * w_c^(1-s) * exp(-w/w_c)
    return 2.0 * alpha * (w**s) * (w_c**(1.0 - s)) * np.exp(-w / w_c)

def discretize_bath(N_b, s, alpha, w_c):
    # 采用对数-线性混合密度进行离散化 rho(w) = (N_b+1)/w_c * exp(-w/w_c)
    # 计算积分：j = (N_b+1) * (1 - exp(-w_j/w_c)) => w_j = -w_c * ln(1 - j/(N_b+1))
    frequencies = []
    couplings = []
    
    # 避免边界奇点，j 从 1 到 N_b
    for j in range(1, N_b + 1):
        # 求解频率 w_j
        w_j = -w_c * np.log(1.0 - j / (N_b + 1.0))
        frequencies.append(w_j)
        
        # 局部态密度 rho(w_j)
        rho_wj = ((N_b + 1.0) / w_c) * np.exp(-w_j / w_c)
        
        # 谱密度
        J_wj = get_sub_ohmic_spectral_density(w_j, s, alpha, w_c)
        
        # 耦合强度 c_j
        c_j = np.sqrt((2.0 / np.pi) * w_j * J_wj / rho_wj)
        couplings.append(c_j)
        
    return np.array(frequencies), np.array(couplings)

def run_sbm_ttns_simulation(s, alpha, N_b=1000):
    # 物理参数
    Delta = 1.0
    w_c = 10.0 * Delta
    t_max = 50.0
    dt = 0.05
    n_steps = int(t_max / dt)
    
    # 获取离散化的浴模式参数
    freqs, couplings = discretize_bath(N_b, s, alpha, w_c)
    
    # 计算绝热重整化隧穿元 Delta_eff (由于 N_b=1000 已经很大，此修正项非常精细)
    # 本处简化演示直接使用 Delta，在严格复现中需积分修正
    Delta_eff = Delta
    
    # 1. 定义物理基底 (Basis Set)
    # 自旋基底: Half-Spin 两能级系统
    spin_basis = BasisHalfSpin(id="spin")
    
    # 玻色子谐振子基底: 截断物理维度 d = 16
    boson_bases = []
    for idx, w_j in enumerate(freqs):
        boson_bases.append(BasisOscillator(id=f"mode_{idx}", n_basis=16, omega=w_j))
        
    # 2. 构建哈密顿量符号项 (Operator terms)
    terms = []
    # 自旋项: -0.5 * Delta_eff * sigma_x
    terms.append(Op("-0.5 * Delta_eff * sigma_x", "spin", Delta_eff=Delta_eff))
    
    # 浴自能项及耦合项
    for idx in range(N_b):
        w_j = freqs[idx]
        c_j = couplings[idx]
        mode_id = f"mode_{idx}"
        
        # 浴能级项: 0.5 * p^2 + 0.5 * w_j^2 * x^2 => 转化为标准谐振子算符形式 H_bath = w_j * a_dagger * a
        terms.append(Op("omega * a_dagger * a", mode_id, omega=w_j))
        
        # 线性系统-浴耦合项: 0.5 * sigma_z * c_j * x => x = 1/sqrt(2*w_j) * (a + a_dagger)
        # 算符系数 coeff = 0.5 * c_j * (1/sqrt(2*w_j))
        coeff = 0.5 * c_j / np.sqrt(2.0 * w_j)
        terms.append(Op("coeff * sigma_z * (a + a_dagger)", ["spin", mode_id], coeff=coeff))
        
    # 3. 树张量网络拓扑设计 (Hierarchical Tree Topology Design)
    # 为1001个物理节点（1自旋+1000谐振子）构建一棵多叉分层平衡树
    # 树结构可以将自旋作为根节点，并将模式按照频率高低对称挂载
    tree_structure = Ttns.build_balanced_tree(spin_basis, boson_bases, branch_ratio=4)
    
    # 4. 构建树张量网络算符 (TTNO)
    ttno_hamiltonian = Ttno(tree_structure, terms)
    
    # 5. 准备初始态 (Hartree Product State: |Up> x |0, 0, ..., 0>)
    # 自旋初始处于 sigma_z = +1 激发态
    spin_init = np.array([1.0, 0.0]) 
    # 玻色子谐振子初始均处于真空态 |0>
    boson_inits = [np.eye(1, 16)[0] for _ in range(N_b)]
    
    init_state = Ttns.build_product_state(tree_structure, [spin_init] + boson_inits)
    
    # 6. 配置投影分裂 TDVP 算法，设置最大键维数 max_bond_dim = 64
    tdvp_propagator = TdvpPS(init_state, ttno_hamiltonian, max_bond_dim=64, threshold=1e-6)
    
    # 7. 时间演化与纠缠熵提取
    time_points = []
    sigma_z_values = []
    spin_entropy_values = []
    
    for step in range(n_steps):
        # 向前推进一个时间步 dt
        tdvp_propagator.evolve(dt)
        current_time = (step + 1) * dt
        
        # 提取当前的 TTNS 量子态波函数
        current_psi = tdvp_propagator.state
        
        # 计算局域可观测物理量: <sigma_z>
        sigma_z_val = current_psi.expectation_value(Op("sigma_z", "spin"))
        
        # 计算系统-浴两体剪缩密度矩阵
        reduced_density_matrix = current_psi.get_reduced_density_matrix("spin")
        
        # 数值求解减缩矩阵的特征值以计算冯·诺伊曼纠缠熵
        eigenvalues = np.linalg.eigvalsh(reduced_density_matrix)
        # 过滤数值微小零特征值，避免 log 崩溃
        eigenvalues = eigenvalues[eigenvalues > 1e-12]
        spin_entropy = -np.sum(eigenvalues * np.log(eigenvalues))
        
        # 打印并存储数据
        time_points.append(current_time)
        sigma_z_values.append(sigma_z_val)
        spin_entropy_values.append(spin_entropy)
        
        if step % 20 == 0:
            print(f"Time: {current_time:.2f} a.u. | <sigma_z>: {sigma_z_val:.4f} | S_spin: {spin_entropy:.4f}")
            
    return np.array(time_points), np.array(sigma_z_values), np.array(spin_entropy_values)

if __name__ == "__main__":
    # 运行典型相干边界点测试
    t, sz, s_ent = run_sbm_ttns_simulation(s=0.7, alpha=0.3, N_b=1000)

3.3 物理复现技术避坑指南与 Checklist

在自主复现该工作中，读者必须特别注意以下技术要点：

高阶玻色子占有数截断验证：虽然一般情况下 $d=16$ 足够，但在强谱发散区域（如 $s \le 0.3, \alpha \ge 0.5$），由于自旋-环境强极化位移，极低频谐振子的平均声子占有数可轻松突破 $20$。因此必须时刻监控末态玻色基底最大物理维度的概率分布，如出现边缘概率溢出，须将对应的低频谐振子基底截断扩展至 $d=32$ 或更高。
TTNS 分叉结构设计：避免使用单一的线性链（这会使 TTNS 退化为 MPS 并丧失长处）。推荐使用四叉或八叉平衡树拓扑，并将自旋置于树的核心几何位置，低频谐振子挂载在距离自旋最近的第一层分支，高频谐振子延伸至边缘分支。这符合大系统耦合中“慢变量主导全局”的物理直觉。
时间步长 $\delta t$ 的收敛性校验：TDVP-PS 具有无条件稳定性，但时间步长过大会导致非流形投影误差积累。务必同时测试 $\delta t = 0.05$ 和 $\delta t = 0.01$。若极化曲线重合，可采用较大的步长以加速扫描。
奇异值自适应截断门槛（SVD Threshold）：自适应键维数截断门槛必须设置为 $\le 10^{-6}$。如为追求速度而放宽至 $10^{-4}$，会产生累积能量不守恒，导致纠缠熵计算出现人工虚假漂移。

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

4.1 关键参考文献及其在工作中的脉络关系

本工作的研究深度紧密依赖于以下里程碑式文献的奠基作用：

Leggett, A. J., et al. Reviews of Modern Physics 1987, 59, 1-85.
贡献与脉络: 自旋-玻色模型的理论根基，奠定了耗散量子系统物理机制的整体框架。
Otterpohl, F.; Nalbach, P.; Thorwart, M. Physical Review Letters 2022, 129, 120406.
贡献与脉络: 首次提出基于自旋极化局域流将亚欧姆 SBM 划分为相干、非相干及伪相干动力学区，是本项纠缠结构研究的最直接物理对标起点。
Haegeman, J.; Lubich, C.; Oseledets, I.; Vandereycken, B.; Verstraete, F. Physical Review B 2016, 94, 165116.
贡献与脉络: 发展了投影分裂（Projector-Splitting）TDVP 积分方案，为基于树拓扑的高维动力学变分算法扫清了极硬代数积分障碍。
Wang, H.; Thoss, M. Chemical Physics 2010, 370, 78-86.
贡献与脉络: 系统探讨了采用多层多组态随时间变化哈特里方法（ML-MCTDH）求解亚欧姆自旋-玻色模型的物理离散化与绝热重整化方法，为本文的浴离散技术提供了重要物理参考。
Li, W.; Ren, J.; Yang, H.; Wang, H.; Shuai, Z. The Journal of Chemical Physics 2024, 161, 054116.
贡献与脉络: 研发了自动 TTNO 构建图论算法及开源张量网络包 Renormalizer，是本项计算最关键的算法平台支撑。

4.2 本工作局限性与前沿学术探讨

尽管本项研究在物理图景和数值算法上展现出无与伦比的精确性，但在面对现实中更复杂的化学/物理现象时，其依然存在以下不可忽视的局限性，指明了下一步的前沿研究方向：

1. 严格零温约束（$T=0$ Limit）

局限性描述：本工作的所有计算均基于纯态波函数 $\|\Psi(t)\rangle$，对应的热力学环境被假定在绝对零温下（所有玻色子初始处于真空零振动态 $|0\rangle$）。然而，对于化学反应、生物光合作用光致激发能转移以及室温分子晶体中的电荷传输，环境温度带来的热涨落和初始混态效应（Thermal Fluctuations & Initial Mixed States）是不可忽略的主导物理力量。
未来突破方向：若要在有限温度（$T > 0$）下进行复现或扩展，必须将当前的 TTNS 升级为纯化热态表象（Thermofield Double, TFD）或直接对树张量算符密度矩阵（Tree Tensor Operator, TTO）进行超算符耗散演化。这将使虚拟键维数翻倍，计算负载面临指数级跃升。

2. 对称无偏置系统限制（Symmetric $\epsilon = 0$ Limit）

局限性描述：本项研究集中在无能量偏置的情形。在分子电子转移（Marcus Theory）等实际场景中，反应物与产物态之间通常存在显著的自由能差 $\epsilon \neq 0$。非对称偏置会剧烈打破双势阱的对称性，从而在 Bloch 几何中极大地倾斜或压制三维轨迹，导致平稳纠缠熵景观的几何脊线拓扑发生彻底重组。
未来突破方向：迫切需要扫描含有非零能级偏置的整个 $(s, \alpha, \epsilon)$ 三维物理相空间，以观察伪相干钩状轨迹在强非对称势阱下的湮灭规律。

3. 谐振子浴离散化极限（Poincaré Recurrence Limit）

局限性描述：将连续浴离散化为 $N_b=1000$ 个谐振子，其等效于在一个准周期的多维有限盒子中进行演化。虽然本工作通过庞大的模式数目将庞加莱复归时间推迟到了 $t_{\text{recurrence}} > 50$ a.u. 之后，但如果研究者需要探索亚欧姆超慢响应的长时极限行为（例如，探索 $t \to 1000$ a.u. 的极端非马尔可夫演化），离散模式产生的虚假人工干涉反馈波将不可避免地污染物理信号。
未来突破方向：必须在张量树的边缘叶子节点处引入非反射开放边界条件或物理耗散型虚部吸收势，以模拟无限连续环境。

4. 马鞍区（$s \approx 0.4$）的超慢收敛难题

局限性描述：在静态熵景观的马鞍凹陷区（$s \approx 0.4$），自旋纠缠熵表现出极慢的对数漂移，即使在 $t = 50$ a.u. 时依然未能完全收敛。尽管作者通过延长计算到 $t = 200$ a.u. 确认了趋势，但这一区域多体关联的高维混沌特性如何定量解析，仍是遗留的物理难题。

5. 补充深度分析：TTNS 与经典 ML-MCTDH 的方法论对决及化学应用扩展

5.1 树张量网络态（TTNS）与多层多组态随时间变化哈特里方法（ML-MCTDH）的深度对决

在分子动力学和理论化学界，ML-MCTDH 长期以来被视作大自由度非谐振浴动力学模拟的垄断性标杆工具。由于 TTNS 和 ML-MCTDH 在波函数表示上都采用了类似的层级多叉树结构（Hierarchical Tree Structure），学者们经常将两者混淆。在此，我们从底层数学机制与计算性能层面对两者进行深入剖析：

物理维度与算法特征	树张量网络态 + TDVP-PS (本文方法)	经典多层 MCTDH (ML-MCTDH)
底层波函数表象	采用离散的正交分域张量分解表示，节点间通过虚拟物理键（Bonds）传递纠缠关联。	采用时变的正交基函数（SPFs）的层级乘积形式构建高维时空波函数。
变分控制原理	将薛定谔方程严格投影至常数虚拟键维数（或自适应精度）的非线性流形切空间上。	基于狄拉克-弗伦克尔变分原理推导，SPFs 本身随时间自适应演化。
数值奇异性问题	无任何奇异性风险。投影分裂（PS）算符演化方案天然避开了密度矩阵求逆，在极低纠缠或态接近纯态时运行极其稳健。	极高奇异性风险。其运动方程包含单粒子减缩密度矩阵的逆矩阵（$\rho^{-1}$）。当某模式的布居数极低（接近纯态或初始零温演化伊始），$\rho$ 的特征值趋于 $0$，方程直接发生零分母崩溃，被迫引入极其敏感的人工正则化参数。
虚拟维数控制	虚拟物理键维数（Bond Dimension $M$）可以沿每条链接独立进行动态 SVD 截断优化，灵活性极大。	每一层子空间的单粒子物理基矢数目（$N_j$）在传播开始前必须人工静态配好，极难做到长时全自适应。
对极端宽带谱稳定性	极强。通过 PS 算符分裂，将具有悬殊能量阶跃的哈密顿量完全解耦成局部平缓物理演化，步长不受高频谱制约。	较弱。由于存在非线性项耦合，常微分积分器（如 VMF 框架）在遇到高频模式振荡时，步长会骤减，导致传播极慢。

总结： TTNS 凭借其在切空间上的线性投影分裂（TDVP-PS），彻底消除了困扰经典分子动力学半个世纪的“减缩密度矩阵求逆奇异性”噩梦。这使其在处理如零温起始、长时低纠缠向高纠缠自然演化等物理过程时，拥有了无可比拟的代数稳定性与超长时步传播速度。

5.2 亚欧姆自旋-玻色模型在复杂化学体系中的微观映射

本研究揭示的亚欧姆物理机制，对理解凝聚态化学和功能材料内部的多体输运动力学具有极强的现实指导意义。

1. 有机半导体中的极化子输运（Polaron Transport in Organic Semiconductors）

在高度各向异性的有机共轭分子晶体（如红荧烯 Rubrene、并五苯 Pentacene）中，载流子（电子或空穴）通常与晶格的低频声学支声子（Acoustic Phonons）和分子间剪切振动发生极强的非绝热耦合。这些低频集体声子运动的谱指数 $s$ 往往在 $0.1 \sim 0.5$ 范围内，呈现出微观的亚欧姆关联特征。

物理映射：自旋代表电子在近邻分子间的跃迁状态（$|\text{Donor}\rangle$ 与 $|\text{Acceptor}\rangle$），而玻色子集合则对应声学支晶格振动。
纠缠启示：本研究发现的“相干相中存在显著中介模间关联（$S > I$）”表明，如果能够通过分子剪裁将电子输运控制在相干隧穿机制内，电子作为相干中介，将会在空间上相隔较远的分子间声子模式间建立协同多体量子关联。这种间接关联能够重塑晶格形变场（Polaron Cloud），从而从根本上降低极化子输运的有效电荷复合能垒。

2. 生物光合复合物中的高效激子能量传递（Coherent Exciton Transfer in Photosynthetic Complexes）

在 FMO 捕光复合物中，激子在不同细菌叶绿素分子（BChls）间的超快能量传递效率在室温下高达 95% 以上。近年来物理界广泛争论的“室温量子相干性”是否扮演了核心角色，可以在此找到纠缠层面的物理解释。

纠缠启示：FMO 复合物外部蛋白质脚手架（Protein Scaffold）所提供的复杂非平衡振动环境常呈现低频亚欧姆特征。本工作指出的“自旋纠缠熵快速锁定至平稳高原”这一现象，解释了为什么即使分子环境存在极强的非马尔可夫噪声，激子隧穿在长时依然能保持极高纯度（$r$ 保持在适度范围）：因为激子在极早期就已经同周围的蛋白骨架振动纠缠饱和，形成了一个宏观稳定的极化子系统，使其免受之后环境细微热涨落的随机去相位干扰。

6. 结论

本项工作通过树张量网络态结合投影分裂时间步变分（TTN-TDVP-PS），不仅全面地勾勒出了亚欧姆自旋-玻色模型在整个 $(s, \alpha)$ 参数面上的底层量子纠缠景观，而且成功地构建起了一条连接“宏观单体布居数分类”与“微观多体系统-环境量子纠缠网络”的物理桥梁。研究指出的“时空尺度分离现象”、“大 $s$ 分叉缺失”、“Bloch 几何钩状轨迹”以及“自旋中介环境关联”，极大地深化了我们对强退相干物理本质的理解，也彰显了现代张量网络技术在解决化学物理经典多体非平衡动力学难题上的澎湃威力。