量子多体引力模拟中的 Krylov 复杂度研究

这是本科毕设（量子多体模型的 Krylov 复杂度研究）的一个简短版本的 note（其实就是草稿），摘要太长了就不放了，简单来说：量子混沌最近出了个新的量化标准——Krylov 复杂度，以及伴随的普适算符增长假说，另一方面，最近发现量子模拟中弯曲时空中的无质量量子场方程可以直接对应一个量子多体模型——跳跃振幅随时间变化的一个最近邻跳跃模型。于是我们将二者结合，在新的模型中计算了量子混沌的一些量，没了，要说有点意思的创新点的话，在数值计算中做的各种优化是相当累人的，也是最令我满意的，至于这个量子模拟模型配合 Krylov 复杂度到底有什么意义，我不知道啊（

KEY WORDS: Quantum chaos; Quantum many-body model; Quantum simulation; Quantum field theory

量子混沌部分主要参考文献是 “普适算符增长假说”¹，弯曲空间量子场映射到量子多体模型的参考文献是 理论部分² 和 实验部分³

没废话的引言

上面已经说了一部分了，那下面稍稍再说一点吧（

量子混沌不是经典混沌中蝴蝶效应的直接对应，因为蝴蝶效应也就是初值敏感性要求系统是非线性的，量子力学线性可积所以做不到，故这里抛弃波函数去考虑算符动力学，然后考虑小扰动在系统是如何传播的就行。比较常用的量化是 OTOC 即非时序关联函数，然后“普适算符增长假说”¹这篇文章又提出了个新的假说和 Krylov 复杂度（其实还提出了更为广义的 Q-复杂度），嘛，二者谁更好也说不清

至于弯曲空间量子场在量子多体系统中的模拟嘛，我不好说，这篇理论文章²提出了个看起来挺对的映射关系，然后解释了一下 Hawking 辐射什么的，又算了算 OTOC 和内外的纠缠熵，呃……说真模型中出现的隧穿效应就是因为离散模型的转变点无法严格为零啊，导致概率不大但是波函数确实可以透射过去（这倒是挺有意思的，这似乎在说黑洞辐射本质是位置离散导致的量子隧穿？），至于 OTOC 算出来发现黑洞是混沌的，这真不好说啊，和本文用 Krylov / Lanczos 方法算出来的不太一样啊（不知道是不是初始算符取的不一样的原因）

Krylov 复杂度

算符空间与 Lanczos 算法

我们的研究对象是算符的演化，所以采用了 Heisenberg 绘景进行描述，算符 \(\mathcal{O}\) 的演化由 Heisenberg 动力学方程 $\dot{\mathcal{O}} = \ii \comm{H}{\mathcal{O}}$ 描述。下面，我们将算符作为一种矢量，称作算符态，记为 $\oket{\mathcal{O}}$，并定义作用于算符态的算符为超算符，显然， \(\require{mathtools} \mathcal{L} \coloneqq \comm{H}{\cdot}\) 是超算符，或称 Liouvillian 超算符。随着时间的推移，算符态会在 Hilbert 空间中运动，在此，我们赋予无限温度内积为 $\obraket{\mathcal{O}_1}{\mathcal{O}_2} \coloneqq \Tr[\mathcal{O}_1^\dagger \mathcal{O}_2] / \Tr[1]$，同理，模长为 $\norm{\mathcal{O}} \coloneqq \obraket{\mathcal{O}}{\mathcal{O}}^{1/2}$。

为了考查算符态的演化，我们将任意时刻的算符态形式解在时间上展开为级数，有 \begin{equation} \oket{\mathcal{O}(t)} = \ee^{\ii \mathcal{L} t} \oket{\mathcal{O}} = \oket{\mathcal{O}} + \frac{\ii t}{1 !} \oket{\mathcal{L} \mathcal{O}} + \frac{(\ii t)^2}{2 !} \oket{\mathcal{L}^2 \mathcal{O}} + \cdots \text{。} \end{equation} 可以发现，在算符态的全部演化过程中，其 Hilbert 空间由序列 \(\{\oket{ \mathcal{L}^n \mathcal{O}}\} = \{\oket{\mathcal{O}}, \oket{\mathcal{L} \mathcal{O}}, \oket{ \mathcal{L}^2\mathcal{O}}, \cdots \}\) 所张成，随着时间的增长，初始算符 $\oket{\mathcal{O}}$ 逐渐扩散到 $\oket{\mathcal{L}^{n} \mathcal{O}}$ 之类的更高阶的项当中，也就是说算符随时间变得更加复杂，这称之为算符增长。现在，因为算符态一定在上述算符态序列张成的空间中运动，所以我们只需要将这个序列正交归一化就可以得到空间的完备基了。这里采用的是 Lanczos 递推算法，其详细步骤如下：

1. 将初始算符态归一化 $\oket{\mathcal{O}_0} = b_0^{-1} \oket{\mathcal{O}}, \ b_0 = \sqrt{\obraket{\mathcal{O}}{\mathcal{O}}}$；
2. 作为递推开始的第一步，利用初始算符态得到 $\oket{\mathcal{O}_1} = b_1^{-1} \mathcal{L} \oket{\mathcal{O}_0}$，其中 $b_1 = \sqrt{\obraket{\mathcal{O}_0 \mathcal{L}}{\mathcal{L} \mathcal{O}_0}}$，若 $b_1 = 0$ 则停止算法；

3. 接着可以开始迭代递推过程，同样若 $b_n = 0$ 则终止迭代，

   \[\begin{equation} \begin{aligned} \oket{A_n} &= \mathcal{L} \oket{\mathcal{O}_{n-1}} - b_{n-1} \oket{\mathcal{O}_{n-2}} \text{，} \\ b_n &= \sqrt{\obraket{A_n}{A_n}} \text{，} \\ \oket{\mathcal{O}_{n}} &= b_n^{-1} \oket{A_n} \text{。} \end{aligned} \nonumber \end{equation}\]

经过上述算法，我们得到了一个正数序列 \(\{b_n\}\) 和一组正交归一基底 \(\{\oket{\mathcal{O}_n}\}\)，前者称为 Lanczos 系数，后者即 Krylov 基。若算法在某处终止迭代，则算符的 Hilbert 空间是有限维的。这里需要注意到的是 Krylov 基只是相对于初始算符 $\oket{\mathcal{O}}$ 而言是完备的，可以描述其经过任意时间的演化后的算符态，但是对于哈密顿量所在的全 Hilbert 空间不一定是完备的，故 Krylov 基张成的空间也被称为 Krylov 子空间。

当有了基底，我们就可以对 Liouvillian 超算符进行展开，其矩阵形式如下： \(\begin{equation} L_{nm} = (\mathcal{O}_n|\mathcal{L}|\mathcal{O}_m) = \begin{bmatrix} 0 & b_1 & 0 & 0 & \cdots \\ b_1 & 0 & b_2 & 0 & \cdots \\ 0 & b_2 & 0 & b_3 & \cdots \\ 0 & 0 & b_3 & 0 & \ddots \\ \vdots & \vdots & \vdots & \ddots & \ddots \end{bmatrix} \text{。} \end{equation}\) 同样我们可以对于算符态进行展开，有 $\oket{\mathcal{O}(t)} = \sum_n \ii^n \phi_n(t) \oket{\mathcal{O}_n}$，其中系数 $\phi_n(t) = \ii^{-n} \obraket{\mathcal{O}_n}{\mathcal{O}(t)}$。

半无限长链与 Krylov 复杂度

我们所得到的系数 $\phi_n(t)$ 实际上是半无限长链的波函数，这一点可以将定义代入 Heisenberg 方程得出，即得离散的 Schrödinger 方程 \begin{equation} \pp_t \phi_n = - b_{n+1} \phi_{n+1} + b_n \phi_{n-1} , \quad \phi_n(0) = \delta_{n0}, \quad n = 1, 2, 3, \cdots \text{。} \label{math:1.1} \end{equation} 容易发现这代表了紧束缚最近邻近似下的跳跃问题，其中跳跃振幅依赖于格点位置。此外，波函数满足 $\sum_n \abs{\phi_n(t)}^2 = 1$，即归一化条件。

算符态从初始时刻开始演化，在半无限长链的角度看来就是位于初始位置的粒子向远处跳跃，于是我们对于算符变得复杂这一抽象的概念有了一个形象的图像，那量化这个图像的关键就是粒子在半无限长链上的平均位置，即 \begin{equation} K(t) = \sum_{n=1} n \abs{\phi_n(t)}^2 \text{，} \end{equation} 这称为 Krylov 复杂度。

除了 Krylov 复杂度以外，另一个比较重要的研究对象是自关联函数，它被定义为算符与初始算符的内积，即 \begin{equation} C(t) = \Tr[\mathcal{O}^{\dagger}(0) \mathcal{O}(t)] / \Tr[1] = (\mathcal{O}|\exp(\ii \mathcal{L} t)|\mathcal{O}) \text{，} \end{equation} 而在半无限长链模型和 Krylov 基展开下，它有一个更加好用的形式： \begin{equation} C(t) = (\ee^{\ii L t})_{00} = \phi_0(t) \text{。} \end{equation}

有关量子混沌与 Krylov 复杂度，Daniel 提出了一个普适的假说¹：在量子混沌系统中，Lanczos 系数应该尽可能快地增长，且最大增长速率为线性，另外在一维系统需要添加一个对数修正。如果用数学的语言来描述，那么对于一个维数为 $d$ 的无限、不可积的量子多体系统而言，其 Lanczos 系数的渐近行为如下： \(\begin{equation} b_n = \begin{cases} A \frac{n}{\ln{n}} + \mathrm{o}(\frac{n}{\ln{n}}) , \quad d=1 \\ \alpha n + \gamma + \mathrm{o}(1) , \quad d>1 \end{cases} \text{。} \end{equation}\) 式中，$\alpha>0$、$A>0$ 和 $\gamma$ 是一些实数常数。

本章后面先是数值验证了“普适算符增长假说”（基于一维横纵场 Ising 模型），其实就是复现参考文献中的那个图

然后说明了 Lanczos 迭代算法存在不稳定性，介绍了两个改进算法，完全正交化算法和部分正交化算法，嘛，这其实是我们数值计算时怀疑自己的计算结果不正确，找了半天文献发现此算法其实数值不稳定，非常高兴啊，以为找到自己哪里算错了（事实上来说，前 100 个 Lanczos 系数一般不存在精度问题哦 XD

弯曲时空量子模拟

弯曲时空量子场

我们的目标是将一个用来描述黑洞的弯曲时空度规映射到量子系统当中，故不失一般性，取 1+1 维 Schwarzschild 坐标 $(t, x)$ 下度规线元为 \begin{equation} \dd{s}^2 = f(x) \dd{t}^2 - f(x)^{-1} \dd{x}^2 \text{。} \label{math:2.1} \end{equation} 我们假定黑洞是静态的且只有一个非退化事件视界，因此要求在视界外部 $x > x_h$ 有 $f(x) > 0$，而且只有一个点 $x = x_h$ 满足 $f(x_h) = 0$ 和 $f’(x_h) > 0$。

因为视界的存在，这个时空存在奇异性，但是其本质上是坐标奇异性，我们选取另外的坐标系是可以消除坐标奇点的。为此，做坐标变换 \begin{equation} v = t + \int f(x)^{-1} \dd{x} \text{，} \end{equation} 我们就能转换到超前 Eddington-Finkelstein 坐标 $(v, x)$ 中，这时度规线元 \(\eqref{math:2.1}\) 变为 \(\begin{equation} \dd{s}^2 = f(x) \dd{v}^2 - 2 \dd{v}\dd{x} \text{，} \label{math:2.2} \end{equation}\) 同时，为了之后的计算方便，在此给出度规张量的矩阵系数： \(\begin{equation} g_{\mu \nu} = \mqty[f(x) & -1 \\ -1 & 0] , \quad g^{\mu \nu} = \mqty[0 & -1 \\ -1 & -f(x)] \text{。} \nonumber \end{equation}\)

考虑弯曲时空中的 Klein-Gordon 方程 $\nabla^{\mu} \nabla_{\mu} \phi + m^2 \phi = 0$，对于复标量场 $\phi$ 和度规 \(\eqref{math:2.2}\)，方程可有具体形式： \(\begin{equation} \begin{aligned} 0 &= \nabla^{\mu} \nabla_{\mu} \phi + m^2 \phi \\ &= g^{\mu \nu} \pp_{\nu} \pp_{\mu} \phi + g^{\mu \nu} \Gamma^{\lambda}_{\nu \mu} \pp_{\lambda} \phi + m^2 \phi \\ &= -f(x) \pp_{x}^2 - 2 \pp_{x} \pp_{v} \phi - f'(x) \pp_{x} \phi + m^2 \phi \text{，} \end{aligned} \label{math:2.3} \end{equation}\) 其中，\(\Gamma^{\lambda}_{\nu \mu} = \frac{1}{2} g^{\lambda \sigma} \qty(\pp_{\mu}g_{\sigma \nu} + \pp_{\nu}g_{\mu \sigma} - \pp_{\sigma}g_{\nu \mu})\) 是 Christoffel 记号，\(\nabla_{\mu}\) 代表协变导数。

为了让方程线性化，我们引入一个新的变量 $\varphi$，其通过 $m \varphi = 2 \pp_{v} \phi + f \pp_{x} \phi$ 来定义，在这个变换下，方程 \(\eqref{math:2.3}\) 转变为一组耦合的一阶方程组： \(\begin{equation} \left\{ \begin{aligned} \pdv{\varphi}{x} &= m \phi \\ \pdv{\phi}{v} &= - \frac{f(x)}{2} \pdv{\phi}{x} + \frac{m}{2} \varphi \end{aligned} \right. \text{。} \end{equation}\) 考虑无质量极限 \(m \rightarrow 0\)，上述两个方程将会解耦，且第一个方程可以被忽略。接着，为了能使连续方程离散化后能得到我们所需要的结果，这里再引入一个变量 $w$，定义为 $w \sqrt{f} = \phi$，此时我们得到了最终的连续线性方程： \(\begin{equation} \pdv{w}{v} = - \frac{1}{4} \qty[\pdv{w f(x)}{x} + f(x) \pdv{w}{x}] \text{。} \label{math:2.4} \end{equation}\)

映射到量子多体系统

下一步，我们将离散化整个系统，取空间坐标为格点位置 $x = x_n = nd, \ n \in \mathbb{N}$，其中 $d$ 为晶格常数，随之有函数离散化为 $f_n = f(nd)$ 和 $w_n(v) = w(v, x_n)$。我们用中心差分来代替空间微分，则方程 \eqref{math:2.4} 变为 \begin{equation} \dv{w_n}{v} = - \kappa_{n+1} w_{n+1} + \kappa_{n} w_{n-1} \text{，} \end{equation} 式中系数 \begin{equation} \kappa_n = \frac{f_n + f_{n-1}}{8d} \approx \frac{f\qty[(n-1/2)d]}{4d} \text{。} \end{equation} 为了让方程的形式看起来更加熟悉，我们再做一个变换 $w_n = \tilde{w}_n (-\ii)^n \ee^{-\ii \mu \nu}$，则上述方程变为 \(\begin{equation} \ii \dv{\tilde{w}_n}{v} = - \kappa_{n+1} \tilde{w}_{n+1} - \kappa_{n} \tilde{w}_{n-1} - \mu \tilde{w}_{n} \text{，} \end{equation}\) 这是一个离散的 Schrödinger 方程，这里系数 $\mu$ 是任意常数。这也代表了一个紧束缚的跳跃问题，对于玻色场或费米场，我们可以很容易地给出二次量子化的哈密顿量： \(\begin{equation} \mathcal{H} = \sum_{n} \qty[-\kappa_n \qty(\hat{a}_n^{\dagger} \hat{a}_{n-1} + \hat{a}_{n-1}^\dagger \hat{a}_n) -\mu \hat{a}_{n}^{\dagger} \hat{a}_{n}] \text{。} \label{math:2.5} \end{equation}\) 其中 $\hat{a}_n^{\dagger}$ 和 $\hat{a}_n$ 分别是第 $n$ 个格点上的产生和湮灭算符，对于玻色子应满足相应的对易关系，而对于费米子则应满足相应的反对易关系。

呃，费米子其实应该用 Dirac 方程来推，当然结果上数学形式是一样的，我们是这样想的，Dirac 方程是 Klein-Gordon 方程的开平方，而这里取了无质量极限，那就没什么区别了

数值计算

在这一部分，我们想在数值上对式 \(\eqref{math:2.5}\) 相似的哈密顿量应用 Krylov 复杂度的相关计算方法，为了使系统变得不可积，在原有的哈密顿量上我们引入了相互作用项，则有新哈密顿量 \(\begin{equation} \mathcal{H} = \sum_{n=2}^{N} \qty[-\kappa_n \qty(\hat{c}_n^{\dagger} \hat{c}_{n-1} + \hat{c}_{n-1}^\dagger \hat{c}_n) + g \hat{c}_{n-1}^{\dagger} \hat{c}_{n-1} \hat{c}_n^{\dagger} \hat{c}_{n}] \text{。} \label{math:3.1} \end{equation}\) 此处的 $g$ 是相互作用强度，$\hat{c}_n$ 和 $\hat{c}_n^{\dagger}$ 分别代表了在第 $n$ 个格点上的费米子的湮灭和产生算符，我们设定系统总共拥有 $N$ 个格点。

为方便进行计算，这里采用 Jordan–Wigner 变换将费米子的产生湮灭算符和 Pauli 自旋算符对应起来，变换的具体形式如下：

\[\begin{equation} \begin{aligned} \sigma_j^{+} &= \ee^{-\ii \pi \sum^{j-1}_{k=1} c_k^{\dagger} c_k} c_j^{+} , \\ \sigma_j^{-} &= \ee^{\ii \pi \sum^{j-1}_{k=1} c_k^{\dagger} c_k} c_j , \\ \sigma_j^{z} &= 2 c_j^{\dagger} c_j - I \text{。} \end{aligned} \nonumber \end{equation}\]

注意到这里我们忽略了算符的帽子，并且 $I$ 是单位算符。在应用了上述变换之后，系统的哈密顿量 \(\eqref{math:3.1}\) 将拥有一个不同的形式： \begin{equation} \mathcal{H} = \sum_{n=2}^{N} \qty[-\kappa_n \qty(\sigma_n^{+} \sigma_{n-1}^{-} + \sigma_{n-1}^{+} \sigma_{n}^{-}) + \frac{g}{4} \qty(\sigma_{n-1}^{z} + I) \qty(\sigma_{n}^{z} + I)] \text{。} \label{math:3.2} \end{equation}

在 Krylov 复杂度的数值计算上，我们的系统有 $N$ 个格点，度规 \(\eqref{math:2.1}\) 中的函数 $f(x)$ 将取双曲正切的形式，则跳跃振幅为 \(\begin{equation} \kappa_{n} = \frac{\beta \tanh \qty[\qty(n - n_h + \frac{1}{2})d]}{4d} \text{。} \end{equation}\) 其中，$n_h$ 是事件视界的位置，另外，我们将 $\mathcal{O}(0) = 2 c_{m}^{\dagger} c_{m} - I = \sigma_{m}^{z}$ 作为了 Lanczos 算法的初始算符。

数值方法与优化

迄今为止，我们还没有详尽地说明过如何计算 Lanczos 系数，对于成为目标的哈密顿量 \(\eqref{math:3.2}\)，将一切在矩阵中展开是合适的。自然起见，取 Pauli 算符的矩阵表示： \(\begin{equation} I = \mqty[1 & 0 \\ 0 & 1] , \ \sigma^{x} = \mqty[0 & 1 \\ 1 & 0] , \ \sigma^{y} = \mqty[0 & -\ii \\ \ii & 0] , \ \sigma^{z} = \mqty[1 & 0 \\ 0 & -1] \text{，} \end{equation}\) 另外有定义 $\sigma^{+} = \qty(\sigma^{x} + \ii \sigma^{y}) / 2 , \ \sigma^{-} = \qty(\sigma^{x} - \ii \sigma^{y}) / 2$。于是对于单个格点而言，其算符可展开在二阶复方阵的线性空间当中，记为 $\mathbb{C}^{2 \times 2}$。我们通常采用算符直积来描述多粒子系统，设系统总共有 $N$ 个格点，则系统所在的线性空间为 $N$ 个 $\mathbb{C}^{2 \times 2}$ 空间的直积，即 $\mathbb{C}^{2^{N} \times 2^{N}}$，对于各个格点上的算符也有同样的直积定义： \begin{equation} P = P_{1} \otimes P_{2} \otimes \cdots \otimes P_{N - 1} \otimes P_{N} \text{，} \end{equation} 则第 $i$ 个格点上的 Pauli 算符为 \(\begin{equation} \sigma^{\alpha}_{i} = \underbrace{I \otimes \cdots \otimes I}_{(i-1) \text{个} I} \otimes \ \sigma^{\alpha} \otimes \underbrace{I \otimes \cdots \otimes I}_{(N - i) \text{个} I} , \quad \alpha \in \{x,y,z,+,-\} \text{。} \end{equation}\)

在 Lanczos 算法中，对于无限温度内积的应用仅限于计算归一化系数，即仅仅计算了算符态的自内积，故存在优化的可能。对于线性空间 $\mathbb{C}^{2^{N} \times 2^{N}}$ ，算符态的 Hilbert 空间中无限温度内积定义为 \(\obraket{\mathcal{O}_1}{\mathcal{O}_2} \coloneqq \Tr[\mathcal{O}_1^\dagger \mathcal{O}_2] / \Tr[1] = \Tr[\mathcal{O}_1^\dagger \mathcal{O}_2] / 2^N\)，故计算过程中所用的自内积可化简为： \(\begin{equation} \obraket{A}{A} = \frac{1}{2^N} \Tr[A^{\dagger} A] = \frac{1}{2^N} \sum_{i, j} A^{*}_{ij} A_{ij} \text{。} \end{equation}\) 式中，$A$ 是任意算符态在线性空间 $\mathbb{C}^{2^{N} \times 2^{N}}$ 中的矩阵表示，$A_{ij}$ 代表了具体的矩阵元，$A^{*}_{ij}$ 代表了矩阵元的复共轭。容易发现，算符态的自内积其实就是其矩阵表示中全体元素的模之和除以维数因子 $2^N$，相较于矩阵乘法来说，经过化简后的表达式其算法时间复杂度更低，且更为容易使用并行化技巧进行计算。

我们的计算存在更多特定的针对性优化，在迭代算法开始的阶段，算符态的复杂程度是相当低的，同时，哈密顿量也是简单的，故二者的矩阵表示是极稀疏矩阵。因所需 Lanczos 系数并不多，故完全可以在计算时采用稀疏矩阵进行存储，这将大大缓解内存的占用，但对于需要计算更多 Lanczos 系数的情况，我们发现其计算速度相较于传统的密集矩阵来说较慢。另外，对于我们所需计算的哈密顿量和初始算符来说，若二者的构成中均不存在 $\sigma^{y}_{i}$ 算符，则其矩阵表示一定是实的，且在往后的计算中也一直是实矩阵，可以用线性空间 $\mathbb{R}^{2^{N} \times 2^{N}}$ 来完全地替代计算，这将在数值上减小一半内存占用并成倍地加快计算速度。

当计算得出前 $n$ 个 Lanczos 系数时，我们便可以对关注重点 Krylov 复杂度和自关联函数进行求解计算了，如 \eqref{math:1.1} 所示，可取前 $n$ 个波函数 $\phi_{i}(t)$ 构成微分方程组，做一个小变换 $\varphi_{i}(t) = \phi_{i}(t) (-\ii)^{n}$ 后可写成更为对称的矩阵形式： \(\begin{equation} \dv{\vb*{\varphi}}{t} = \ii L \vb*{\varphi} = \ii \mqty[ 0 & b_1 & 0 & 0 & \cdots \\ b_1 & 0 & b_2 & 0 & \cdots \\ 0 & b_2 & 0 & b_3 & \cdots \\ 0 & 0 & b_3 & 0 & \ddots \\ \vdots & \vdots & \vdots & \ddots & \ddots ] \vb*{\varphi} \text{，} \end{equation}\) 式中，\(\vb*{\varphi} = \qty[\varphi_0, \varphi_1, \cdots , \varphi_{n-1}]^{T}\) 是 $n$ 维列矢量，$L$ 是 Liouvillian 超算符在 Krylov 基下展开系数的 $(n \times n)$ 矩阵。微分方程的形式解我们能很容易地给出 \(\vb*{\varphi}(t) = \exp(\ii L t) \vb*{\varphi}(0)\)，那问题的本质就归结于矩阵 $L$ 的对角化，这是个实对称矩阵，故存在矩阵 $S$ 使得 $L = S D S^{-1}$，其中 $D$ 是对角矩阵，对角线上是 $L$ 的全部特征值，即 $D = \mathrm{diag}(\lambda_{1}, \lambda_{2}, \cdots, \lambda_{n})$。有初始条件 \(\vb*{\varphi}(0) = \qty[1, 0, \cdots , 0]^{T}\)，现在我们可以给出微分方程的解： \(\begin{equation} \vb*{\varphi}(t) = \ee^{\ii L t} \vb*{\varphi}(0) = S \ee^{\ii D t} S^{-1} \vb*{\varphi}(0) = S \mqty[\dmat{\ee^{\ii \lambda_{1} t}, \ee^{\ii \lambda_{2} t}, \ddots, \ee^{\ii \lambda_{n} t}}] S^{-1} \vb*{\varphi}(0) \text{。} \end{equation}\) 注意到 Krylov 复杂度 $K(t) = \sum_{i=1}^{n-1} i \abs{\varphi_{i}(t)}^2$，而自关联函数 $C(t) = \varphi_{0}(t)$，一切似乎都很美好，但是，我们所计算得出的 Lanczos 系数的有限性限制了演化的正确性，在一段时间以后，有限截断会导致 Krylov 复杂度和自关联函数失真，这是很容易预料到的。

对于自关联函数，我们可以直接对哈密顿量进行严格对角化操作来求解，考虑到其定义为 $C(t) = \Tr[\mathcal{O}^{\dagger}(0) \mathcal{O}(t)] / \Tr[1]$，算符的演化遵从 Heisenberg 运动方程，故有 $\mathcal{O}(t) = U^{\dagger}(t) \mathcal{O}(0) U(t)$，其中 $U(t) = \exp(-\ii H t)$ 为时间演化算符，$H$ 是哈密顿量的表示矩阵。我们研究的量子系统的哈密顿量必然是厄米的，则一定可以进行对角化，故存在矩阵 $P$ 使得 $H = P M P^{-1}$，其中 $M = \mathrm{diag}(\varepsilon_{1}, \varepsilon_{2}, \cdots, \varepsilon_{n})$ 为对角矩阵，对角线上是哈密顿量矩阵的特征值。现在，自关联函数可以表示为 \(\begin{equation} \begin{aligned} C(t) &= \frac{1}{2^N} \Tr[\mathcal{O}^{\dagger}(0) P \ee^{\ii M t} P^{-1} \mathcal{O}(0) P \ee^{-\ii M t} P^{-1}] \\ &= \frac{1}{2^N} \Tr[P^{-1} \mathcal{O}^{\dagger}(0) P \ee^{\ii M t} P^{-1} \mathcal{O}(0) P \ee^{-\ii M t}] \\ &= \frac{1}{2^N} \Tr[Y^{\dagger} \ee^{\ii M t} Y \ee^{-\ii M t}] \text{，} \end{aligned} \end{equation}\) 式中定义了新的矩阵 $Y \equiv P^{-1} \mathcal{O}(0) P$ 进行化简，同时不妨将对角矩阵的矩阵元重新定义，即写为 $\exp(-\ii M t) = \mathrm{diag}(\exp(- \ii \varepsilon_{1} t), \exp(- \ii \varepsilon_{2} t), \cdots, \exp(- \ii \varepsilon_{n} t)) \equiv d = \mathrm{diag}(d_{1}, d_{2}, \cdots, d_{n})$。到此为止，我们仍然需要在每一个时间步长上进行多次矩阵乘法运算和求迹运算，这是非常浪费时间的，实际上对于上述表达式仍可以继续化简，将矩阵乘法拆成分量形式进行推导，我们可以得到： \(\begin{equation} \begin{aligned} \Tr[Y^{\dagger} \ee^{\ii M t} Y \ee^{-\ii M t}] &= \sum_{i,j,k,p,q} Y^{\dagger}_{qp} d^{*}_{pi} Y_{ij} d_{jk} \delta_{kq} = \sum_{i,j,k,p,q} Y^{\dagger}_{qp} d^{*}_{pi} \delta_{pi} Y_{ij} d_{jk} \delta_{jk} \delta_{kq} \\ &= \sum_{i,j,q} Y^{\dagger}_{qi} d^{*}_{i} Y_{ij} d_{j} \delta_{jq} = \sum_{i,j} d^{*}_{i} d_{j} Y^{\dagger}_{ji} Y_{ij} \\ &= \sum_{i,j} d^{*}_{i} d_{j} Y^{*}_{ij} Y_{ij} = \sum_{i,j} d^{*}_{i} d_{j} \abs{Y_{ij}}^2 \\ &= \mqty[d^{*}_{1} & d^{*}_{2} & \cdots & d^{*}_{n}] Y \mqty[d_{1} \\ d_{2} \\ \vdots \\ d_{n}] \text{。} \end{aligned} \end{equation}\) 定义 \(\vb*{d} = \mqty[d_{1} & d_{2} & \cdots & d_{n}]^{T}\)，则最终自关联函数为 \(C(t) = \vb*{d}^{\dagger} Y \vb*{d} / 2^N\)，是一个二次型，注意到式中矩阵 $Y$ 是不含时的，仅仅与初始算符和哈密顿量有关，则在每一个时间步长的计算中，我们仅需要计算一次二次型即可，这是十分快速省时的。

Pauli 串算法

一维自旋链是十分常见的物理模型，数学上也对此颇有研究，在量子信息领域，格点 Pauli 自旋矩阵属于 Clifford 群，记为 \(\mathrm{Cl}_{n}\)，其任意群元的幺正变换仍然是群元，我们的目标在于找到更为方便的数学表示，以便于更好地进行数值计算。首先，让我们用新的符号来定义 Pauli 矩阵： \(\begin{equation} \tau_{00} = \mqty[1 & 0 \\ 0 & 1] , \ \tau_{01} = \sigma^{x} = \mqty[0 & 1 \\ 1 & 0] , \ \tau_{10} = \sigma^{z} = \mqty[1 & 0 \\ 0 & -1] , \ \tau_{11} = \ii \sigma^{y} = \mqty[0 & 1 \\ -1 & 0] \text{，} \end{equation}\) 这种表示下所有的矩阵元皆为实数，至于虚数的部分，则被我们拆到了别处，另外，注意到下标是极其独特的，这是二进制的表示方式，或者说是整数域 $\mathbb{Z}_{2}$。对于格点总数为 $n$ 的系统，格点上的 Pauli 算符，即 $n$ 个 Pauli 算符的张量积，可以被表示为 \(\begin{equation} \tau_{a} = \tau_{v_1 w_1} \otimes \tau_{v_2 w_2} \otimes \cdots \otimes \tau_{v_n w_n}, \quad v,w \in \mathbb{Z}^{n}_{2}, \ a = \mqty[v \\ w] \text{。} \end{equation}\)

定义多粒子 Pauli 串群 $G$，其任意群元的有如下表示： \(\begin{equation} g = \ii^{\delta} (-1)^{\epsilon} \tau_{a} \in G, \quad \delta, \epsilon \in \mathbb{Z}_{2}, \ a \in \mathbb{Z}_{2}^{2n} \text{。} \end{equation}\) 在式中，之所以既有虚数部分也有负数符号部分，就是因为我们想让其指数是容易被二进制表示的，下面将看到，这也会简化数值计算。对于群 $G$ 的定义来说，需要同时定义群乘法才能完整地构成一个群，故我们将算符的表示矩阵乘法作为群乘法，对于群元 $g_1, g_2 \in G$，有如下群乘法运算规则： \(\begin{equation} \begin{aligned} & g_1 g_2 = \ii^{\delta_1} (-1)^{\epsilon_1} \tau_{a_1} \ii^{\delta_2} (-1)^{\epsilon_2} \tau_{a_2} = \ii^{\delta_{12}} (-1)^{\epsilon_{12}} \tau_{a_{12}} , \\ & \delta_{12} \equiv \delta_{1} + \delta_{2} \ (\mathrm{mod}\ 2), \\ & \epsilon_{12} \equiv \epsilon_{1} + \epsilon_{2} + \delta_{1} \delta_{2} + \vb*{v}_{1} \vdot \vb*{w}_{2} \ (\mathrm{mod}\ 2), \\ & a_{12} \equiv a_{1} + a_{2} \ (\mathrm{mod}\ 2) \text{，} \end{aligned} \end{equation}\) 式中，\(a_{1} = \mqty[\vb*{v_{1}} & \vb*{w_{1}}]^{T}, \ a_{2} = \mqty[\vb*{v_{2}} & \vb*{w_{2}}]^{T}\)，上述运算中所有数均遵循有限域 $\mathbb{Z}_{2}$ 的运算，即保证所有数只有 0 和 1 两种取值。

现在，我们对于群 $G = {g_{\alpha}}$ 已经有了完整的定义了，对于有限的粒子数 $n$，有限群的阶数为 $2^{2n + 2}$，但是，我们所想要计算的哈密顿量和算符通常不是单个的群元，而是群元的线性组合，故需要建立群空间和群代数的概念。

此处是群论概念，略，简单来说就是群空间是群元作为基矢的空间，群空间里的向量做运算，详见 群论笔记

下面，我们将说明具体数值方案。

对于二进制计算机来说，有限域 \(\mathbb{Z}_{2}\) 的相关运算是令人欣喜的，可以很方便地将其对应到位运算上。假设 \(a, b, c \in \mathbb{Z}_{2}\)，域加法 \(c \equiv a + b \ (\mathrm{mod}\ 2)\) 直接对应于异或运算 \(c = a \oplus b\)，可简单解释为相同为 0 不同为 1；域乘法 \(c \equiv a b \ (\mathrm{mod}\ 2)\) 直接对应于和运算 $c = a \land b$，可简单解释为全 1 则 1 否则为 0。对于任一群元 $g = \ii^{\delta} (-1)^{\epsilon} \tau_{a} \in G$，注意到其存在 $(2n + 2)$ 个二进制位，包括虚数指标 $\delta$，负号指标 $\epsilon$ 和 $2n$ 个二进制位组成的 \(a \in \mathbb{Z}_{2}^{2n}\)。但考虑到我们使用的是群代数 $R_{G}$ 来进行具体运算，其元素是 $G$ 的群元的线性组合，各个系数是复数，故完全可以忽略虚数指标 $\delta$ 和负号指标 $\epsilon$，仅仅考虑 $a$ 即可。进一步，\(a = \mqty[\vb*{v} & \vb*{w}]^{T}\) 可以看作是两个 $n$ 位二进制数串所编码的，换而言之，将二进制转为十进制后，可以使用两个整数 $v$ 和 $w$ 来表示群空间的所有基。在这里我们来举个例子说明这一点：假设系统有 $n=3$ 个粒子，有算符 \(\mathcal{O} = \sigma^{x}_{1} + \sigma^{z}_{3}\)，注意到 \(\sigma^{x}_{1} = \tau_{01} \otimes \tau_{00} \otimes \tau_{00}\) 以及 \(\sigma^{z}_{3} = \tau_{00} \otimes \tau_{00} \otimes \tau_{10}\) ，则可将其映射到群空间中的元素 \(x = \tau_{(000)_{2}, (100)_{2}} + \tau_{(001)_{2}, (000)_{2}} = \tau_{0, 4} + \tau_{1, 0}\)，式中的标记 $(a)_{2}$ 代表 $a$ 是二进制。

在具体计算的过程中，我们发现对于 Lanczos 算法而言，算符态会越来越复杂，若哈密顿量涉及所有格点位置的 Pauli 算符，则算符态映射到群空间中的元素 $x = \sum_{i} c_i \tau_{v_i w_i}$ 时，其群元几乎在单一指标上遍历所有可能值，即 $v_i$ 和 $w_i$ 能单独完整地取到 $0$ 到 $(2^n - 1)$ 的所有值。也就是说，对于此问题，我们可以采用数组来标记某一指标，不妨定为 $v$，数组的下标就是其取值，而对于另一指标 $w$ 来说，当 $v$ 固定时，$w$ 的可能取值较少，我们更倾向于使用哈希散列映射表来存储以便于快速地取值和合并值，其键为 $w$ 而值为复数系数 $c_i$。

总结来说，在此算法中，我们使用数组嵌套哈希表来存储任何算符，数组下标 $v$ 和哈希表的键 $w$ 可以一一对应群空间的基，即部分群元 $\tau_{v w} $，并且所有算符间运算遵循上述所说的群代数运算，则我们将 Lanczos 算法全部落实到了具体的数据结构上了。在实际计算当中，我们发现此算法实现在粒子数 $n$ 较大或者哈密顿量较为简单时比矩阵算法实现消耗了更少的内存，可以相当快速地计算出前几个 Lanczos 系数，具有不错的应用价值，若系统存在某些对称性可以进一步化简，则在数值计算上将具有更大的优势。但是，在大量迭代之后，由于算符态已经变得足够复杂，Krylov 子空间已经趋于饱和，此算法将在计算速度上处于劣势，且因为哈希表结构的不稳定将导致内存占用波动较大，若需要计算很多 Lanczos 系数，或者在 $n$ 较小的系统上进行计算，则不应该使用此算法。

对对，发现了吗，这就是之前写的一篇文章 Hash Table 性能优化 的理论部分，所谓的 Pauli String 算法

数值计算结果

一开始，我们使用了上述的 Pauli 串算法进行问题的数值求解，取了 20 个格点并且事件视界在几乎中间的位置上，然后在超算上算啊算啊算……心惊胆战，内存不知道什么时候就会炸（哈希表嘛，内存占用真的是玄学

然后发现在前 12 个左右 Lanczos 系数和矩阵算法（13 个格点，事件视界在正中间）是一样的，后面才开始偏离，那不管了，能精确对角化才是好的，因为要得到准确的自关联函数啊（如果用有限截断的 Lanczos 系数来算，那时间演化上必须很短，大约 1~2 个时间单位就飞掉了

Krylov 复杂度就算了吧，取个 100 多个 Lanczos 系数才能算个大概，当然这不重要，因为不管怎么画前期都是指数增长的（

• 从 Lanczos 系数上来看，早期都有点线性增长，后期，如果没有相互作用的话，系数会趋于平缓，甚至下降；有相互作用就会线性增长，而且增长斜率和相互作用强度正相关
• 和实验³可以对上的两点
  • 从自关联函数上来看，开放边界条件（就是闭边界条件）的作用和事件视界一样，都会导致图像上出现反弹峰，很好理解嘛，连续模型上那点系数为零，离散的多体模型呢，周围跳跃振幅很小，所以很难传播，一点透射加大量反弹呗
  • 初始算符在事件视界上时，自关联函数下降特别缓慢，这当然也很好理解，跳跃振幅小所以很难跳出去嘛
• 相互作用强度很大时，事件视界的作用几乎消失了，自关联函数看不到反弹峰了，这两者说不定有竞争关系，当然相互作用强度不大时二者皆可存在，所以……这说不定可以在弯曲度规对应的量子多体模型中引入小小的相互作用（即让系统混沌）来研究，新方向？

上面就是结论，各种计算结果图就不放了

论文代码详见 本科毕设论文代码