跳转至

连续空间的世界线 Worm 算法

上一节当中,我们将世界线构型的时间连续化,即虚时方向世界面的数目几近于 \(10^8\) 的量级;而世界面的空间坐标是离散晶格,从而能够模拟晶格上的相互作用模型。 然而,对于许多真实体系的问题,例如计算 \(^4{\rm He}\) 气体的超流转变温度时,空间必然不能取得过于离散化。因此,我们可以放宽对时间连续化的程度,即取到尽可能离散化而不失去表达物理现象的能力,通常最高可到 \(10^4\) 的量级;同时对世界面的空间尽可能取连续化坐标。

相比于晶格上的哈密顿量,连续空间的哈密顿量仍然可以分为对角项 \(\hat{H}_0\) 和非对角项 \(\hat{V}\) :

\[ \hat{H} = \hat{H}_0 + \hat{V} \]
\[ \hat{H_0} = \sum_{i<j} U(\vec{r}_{i,j}), \quad \hat{V} = - \sum_i \lambda_i \nabla_i^2 \]

其中非对角项 \(\hat{V}\) 当中 \(\lambda_i = {\hbar^2}/{2m_i}\) 与不同粒子的质量相关,描述了动能项。对角项 \(\hat{H}_0\) 则描述了不同粒子之间的相互作用。

配分函数空间 \(\mathcal{Z}\)

同样地,我们可以写出连续空间体系的配分函数:

\[ \mathcal{Z} = \sum_{\lbrace \alpha^{(i)} \rbrace} \prod_{i=0}^{N-1} \langle \alpha^{(i+1)} | e^{-\tau \hat{H}} | \alpha^{(i)} \rangle =\sum_{\lbrace \alpha^{(i)} \rbrace} W(\lbrace \alpha^{(i)} \rbrace) \]

其中基矢量 \(\alpha\) 由各个粒子的空间位置组成, 对配分函数其中的一个重复单元,利用 Trotter 分解有:

\[ \langle \alpha^{(i+1)} |e^{-\tau \hat{H}} | \alpha^{(i)} \rangle \approx e^{-\tau H_0} \langle \alpha^{(i+1)} | e^{-\tau \hat{V}}| \alpha^{(i)}\rangle \]

注意这里与之前不同,我们可以解析地写出动能的权重。利用矩阵的直积公式:

\[ (\textbf{A}_1 \otimes \textbf{B}_1) (\textbf{A}_2 \otimes \textbf{B}_2) = (\textbf{A}_1 \textbf{A}_2) \otimes (\textbf{B}_1 \textbf{B}_2) \]

再根据定义 \(| \alpha^{(i)}\rangle = \bigotimes_{j}^{N_p} | \alpha^{(i)}_j\rangle\),将多个粒子的动能拆成单个粒子的乘积:

\[ \langle \alpha^{(i+1)} | e^{-\tau \hat{V}}| \alpha^{(i)}\rangle = \prod_j^{N_p} \langle \alpha^{(i+1)}_j | e^{-\tau \hat{V}}| \alpha^{(i)}_j\rangle \]

其中 \(N_p\) 表示总的粒子数,下指标 \(j\) 表示不同的粒子,上指标 \((i)\) 表示不同的世界面。这里我们相当于将基矢量 \(| \alpha^{(i)}\rangle\) 拆成了不同粒子的直积,即 \(| \alpha^{(i)}_j\rangle\) 的直积。值得注意的是,对于单个粒子的动能权重,\(\hat{V}\) 并不会像之前晶格体系当中那样会筛选掉不满足的态,也就是说不会对下一个世界面当中粒子的位置做出任何限制,权重均由两世界面上粒子位置 \(\vec{r}_j^{(i)}\)\(\vec{r}_j^{(i+1)}\) 来决定,类似于自由粒子的传播子,如下式所示:

\[ \langle \alpha^{(i+1)}_j | e^{-\tau \hat{V}}| \alpha^{(i)}_j\rangle = \left(\frac{1}{4 \pi \lambda_j \tau} \right)^{d/2} \exp{\left( - \frac{(\vec{r}_j^{(i+1)} - \vec{r}_j^{(i)})^2}{4\lambda_j \tau} \right)} \]

其中 \(d\) 表示世界面的空间维度。上式也称为密度矩阵,可以用 \(\rho_0(\vec{r}_j^{(i)},\vec{r}_j^{(i+1)},\tau)\) 表示。至此,我们可以写出配分函数空间权重的形式为:

\[ W(\lbrace \alpha^{(i)} \rbrace) \approx \sum_{\lbrace \alpha^{(i)} \rbrace} \prod_{i=0}^{N-1} \langle \alpha^{(i+1)} | e^{-\tau \hat{V}}| \alpha^{(i)}\rangle \exp \left({-\sum_{i=0}^{N-1} H_0^{(i)} \tau} \right) \]

注意,由于空间是近乎连续的,每个世界面的粒子构型几乎都与下一个世界面不同,因此不需要像之前晶格上那样以 Kink 的数目 \(\mathcal{K}\) 进行分类,而是非对角项和对角项一样,遍历整个虚时。

格林函数空间 \(\mathcal{G}\)

类似地,为了进行 Worm 的更新我们这里可以引入产生湮灭的场算符,定义松原格林函数:

\[ G(I,M,\tau_I,\tau_M) = \left \langle \hat{b}_I(\tau_I) \hat{b}^\dagger_M(\tau_M) \right \rangle = \frac{g(I,M,\tau_I,\tau_M)}{\mathcal{Z}} \]

其中 $ g(I,M,\tau_I,\tau_M) =
{\rm Tr}[\mathcal{T} \hat{b}_I(\tau_I)\hat{b}^\dagger_M(\tau_M) e^{-\beta\hat{H}}] $. 同样地,我们也可以定义格林函数空间的配分函数 \(\mathcal{G}\), 从而更新构型的算法可以在其中抽样:

\[ \mathcal{Z}_w = \mathcal{Z} + \omega_G \mathcal{G} \]

<!--

\[ \mathcal{G} = \sum_{p_I,p_M = 0}^{N-1} \sum_{I,M=0}^{N_s-1} g(I,M,p_I\tau,p_M\tau) \tau^2 = \sum_{p_I,p_M = 0}^{N-1} \sum_{I,M=0}^{N_s-1} W_{\mathcal{G}} \]

-->

注意,由于连续空间当中,粒子严格处于同一个位置的概率几乎为零,因此同一个位置的粒子数最大到是 1,进而因为算符作用产生额外的系数项 \(\sqrt{n_I(n_M+1)} = 1\),与硬核玻色子类似,故可以不用考虑。

因此权重 \(W_{\mathcal{G}}\)\(W\) 的差别仅仅是多了求和变量 \(\tau^2 d\vec{r}_I d\vec{r}_M\),以及要扣除被湮灭掉的或加上被产生出的自由粒子的传播子 \(\rho_0(\vec{r}_j^{(i)},\vec{r}_j^{(i+1)},\tau)\) 所占的权重。

注意,这里的权重 \(W_{\mathcal{G}}\) 与连续时间算法的情况有微妙的不同。连续时间算法当中非对角项 \(\hat{V}\)\(\hat{b}_I\) (或 \(\hat{b}_M^\dagger\))几乎不可能处于同一时空,一方面是因为虚时间是准连续的,出现的概率趋于零;另一方面,在算法的更新上会规避掉同时出现的情况。然而,连续空间算法则不同,这源于其几乎每一个世界面当中的每一个粒子在一下虚时的世界面当中都会有变动,从而虚时上几乎处处都有 \(\hat{V}\) 的贡献, 因此 \(\hat{b}_I\) (或 \(\hat{b}_M^\dagger\))就无法避免地与 \(\hat{V}\) 共存,其结果就是 \(\sqrt{n_I(n_M+1)}\)\(\rho_0(\vec{r}_j^{(i)},\vec{r}_j^{(i+1)},\tau)\) 绑定,因此格林函数空间的权重需要增减 \(\rho_0(\vec{r}_j^{(i)},\vec{r}_j^{(i+1)},\tau)\) 项。

常见体系

\[ \hat{H} = -\frac{1}{2}\sum_{i=1}^N \nabla_i^2 + \frac{1}{2} \sum_{i=1}^N \frac{z_i^2}{\xi^4} + \sum_{i<j} \left\lbrace V_{\rm{dd}}(r_{ij}) + \left(\frac{\sigma}{r_{ij}}\right)^{12} \right\rbrace \]

构型的更新

连续空间 Worm 算法的更新仍然在

\[ \mathcal{Z}_w=\mathcal{Z}+\omega_G\mathcal{G} \]

中进行。与连续时间格点模型不同,这里每个世界面上的粒子坐标都是连续变量,相邻世界面之间的传播子

\[ \rho_0(\vec r,\vec r',\tau) \]

本身就是权重的重要组成部分。因此更新不只是移动虚时位置,还要改变粒子坐标、开放或闭合世界线、增减粒子段,并处理交换置换。

下面只说明各更新的作用和接受率结构。具体实现中,接受率统一由

\[ P_{\rm acc}(\nu\to\nu') = \min\left[ 1, \frac{W(\nu')A(\nu'\to\nu)} {W(\nu)A(\nu\to\nu')} \right] \]

给出;如果构型处于格林函数空间,还需要包含 \(\omega_G\)、worm 端点积分测度和相应的 proposal 概率。

Wiggle/Displace

WiggleDisplace 都在配分函数空间 \(\mathcal Z\) 中更新闭合世界线,不改变粒子数,也不引入 worm 端点。

Wiggle 通常选择某条世界线的一段连续虚时区间,固定两端坐标,只重新抽样中间若干世界面的粒子坐标。它的目的类似于连续时间算法中的局域片段移动:在不改变拓扑结构的情况下,让路径在空间中充分松弛。由于动能传播子是高斯形式,proposal 常参考自由粒子 Brownian bridge,使动能权重尽量被 proposal 吸收,接受率主要由相互作用能变化决定。

Displace 则把整条世界线或一段路径整体平移:

\[ \vec r_i^{(k)}\to \vec r_i^{(k)}+\Delta\vec r. \]

它改变的是粒子在盒子中的整体位置,而不是路径形状。若采用周期边界,需要把坐标重新映射回模拟盒。该更新对液体、气体和弱束缚体系很有用,可以避免世界线质心移动过慢。

这两类更新对应连续时间算法中“移动已有片段”的思想,但连续空间中移动的是坐标轨迹,而不是晶格上的占据数片段。

Open/Close

OpenClose 负责在配分函数空间 \(\mathcal Z\) 与格林函数空间 \(\mathcal G\) 之间切换。

Open 从一条闭合世界线中选取一个虚时片段,将其切开,形成两个 worm 端点。切开后,中间被移除或改写的传播子不再作为闭合路径的一部分,而是进入格林函数空间的权重。两个端点可以理解为连续空间版本的产生、湮灭算符位置。

Close 是反向操作:当两个 worm 端点属于可闭合的同一条开放世界线时,尝试重新连接它们,恢复闭合路径并回到 \(\mathcal Z\) 空间。

接受率中需要比较:

  • 闭合路径权重与开放路径权重;
  • 被删除或补回的自由传播子 \(\rho_0\)
  • 相互作用能积分的变化;
  • 选择世界线、选择虚时区间、选择新坐标路径的 proposal 概率;
  • \(\mathcal G\) 空间相对权重 \(\omega_G\)

在结构上,它对应连续时间算法里的 Create/Delete Worm,但连续空间中还要处理端点之间路径坐标的连续抽样。

Insert/Remove

InsertRemove 用于在格林函数空间中增减一条完整开放世界线段,通常与 grand-canonical ensemble 中粒子数涨落有关。

Insert 尝试在某个虚时区间插入一段新的开放路径。proposal 需要选择起止虚时、端点坐标,并根据自由传播子或 Brownian bridge 生成中间坐标。新路径带来动能传播子、相互作用能和化学势项的变化。

Remove 是反向操作:选择已有开放路径段并删除。若删除后构型满足目标空间的约束,则按反向权重比接受。

这类更新的作用类似连续时间算法中创建或删除一段 segment,但连续空间中 segment 不再只是占据数常数片段,而是一条带有空间坐标轨迹的路径。

接受率必须包含插入路径的积分测度。实现时最容易出错的地方,是 proposal 中坐标采样密度与权重中的自由传播子是否一致;若两者设计得好,许多高斯传播子因子可以在比值中抵消。

Advance/Recede

AdvanceRecede 在格林函数空间 \(\mathcal G\) 中移动 worm 端点的虚时位置,从而改变开放世界线的长度。

Advance 让某个 worm 端点沿虚时方向前进一段距离,并为新增长的路径段生成空间坐标。新增部分的权重来自自由传播子、相互作用能以及可能的化学势项。

Recede 则让端点后退,删除一段已有路径。它是 Advance 的反向过程。

这对应连续时间算法中 worm head 在虚时方向的 Shift Worm,但连续空间中端点移动不仅改变时间长度,也引入或删除一段连续坐标路径。因此 proposal 通常需要同时抽样时间长度和空间轨迹。

接受率中应检查:

  • 新增或删除路径段的传播子权重;
  • 对角相互作用能积分的变化;
  • 端点时间长度 proposal 的正反向概率;
  • 空间路径 proposal 的正反向概率;
  • 周期虚时边界下是否跨过 \(\beta\)

如果端点移动跨过其他粒子世界线或周期边界,构型的置换结构也可能改变,需要与 Swap 更新配合。

Swap

Swap 用于处理相同粒子的置换。对于玻色体系,配分函数需要对所有粒子置换求和;在路径积分表示中,这对应世界线可以在虚时周期边界处互相连接,形成交换环。

局域坐标更新通常很难改变置换拓扑。Swap 的作用是让 worm 端点尝试接到另一条世界线上,从而改变开放路径与闭合路径的连接关系。若最终闭合,构型可能形成不同长度的 permutation cycle。

可以把它理解为连续空间版本的“worm 在空间方向移动”:连续时间格点 Worm 中,worm head 通过创建或删除 kink 从一个格点移动到相邻格点;连续空间中没有固定格点,worm 端点通过选择附近粒子的世界线并重连路径来改变粒子标签和置换结构。

接受率需要比较重连前后的传播子、相互作用能、选择目标世界线的概率,以及生成连接路径的 proposal 概率。对玻色体系,长交换环与超流响应密切相关,因此 Swap 是否高效会直接影响低温超流模拟的采样质量。

补充

  • 首先在配分函数空间当中进行预模拟,一方面根据预模拟结果确定各个操作的接受概率以及level,化学势(不小于焓);
  • 在格林函数空间进行试模拟,保证各个参数指标符合要求,例如pcycle里面尝试次数在 \(10^3\) 量级;
  • 进行正式模拟,可以先对一个温度 \(T\) 下成倍增加虚时层数 \(N\), 确定好哪个 \(N\) 下能量在误差棒下收敛,后续的模拟只需保证 \(N T\) 乘积不变即可;
  • 当温度之外的其他参数发生变化了,以上步骤需要重新进行。