Flowers for Tuesday

Computational Physics

路径积分蒙特卡洛方法 / Path Integral Monte Carlo

下载 Markdown

参考教材

[1] Jos Thijssen , Computational Physics

[2] J.J.Sakurai , Jim Napolitano , Modern Quantum Mechanics

路径积分蒙特卡洛这个名字其实很容易给初学者误导，总让人以为是借助蒙特卡洛方法来解决波函数的路径积分问题。但实际上PIMC的出发点与其关系不大，只是借用了路径积分的方法。

1. 路径积分

不过总还是得先介绍一下什么是路径积分。

路径积分方法是在解决一个这样的问题，假设时刻 tit_i 的初始波函数为 ψ(xi,ti)\psi(x_i, t_i)，我们希望直接计算出任意未来时刻 tft_f 的波函数 ψ(xf,tf)\psi(x_f, t_f)，我们可以定义一个传播子K(xf,tf;xi,ti)K(x_f, t_f; x_i, t_i)满足

ψ(xf,tf)=∫−∞∞K(xf,tf;xi,ti) ψ(xi,ti) dxi\psi(x_f, t_f) = \int_{-\infty}^{\infty} K(x_f, t_f; x_i, t_i) \, \psi(x_i, t_i) \, dx_i

于是等价的，其形式可以写为

K(xf,tf;xi,ti)=⟨xf∣e−iℏH^(tf−ti)∣xi⟩K(x_f, t_f; x_i, t_i) = \langle x_f \vert{} e^{-\frac{i}{\hbar} \hat{H} (t_f - t_i)} \vert{} x_i \rangle

那么现在问题就是如何 将其写成一个完整的数学形式，我们知道哈密顿量为

H^=T^+V^=p^22m+V(x^)\hat{H} = \hat{T} + \hat{V} = \frac{\hat{p}^2}{2m} + V(\hat{x})

注意到动能和势能算符不可对易

[T^,V^]≠0[\hat{T}, \hat{V}] \neq 0

对于非对易算符，eA^+B^≠eA^eB^e^{\hat{A} + \hat{B}} \neq e^{\hat{A}} e^{\hat{B}}。因此这里不能直接将其拆分，而是需要将总时间间隔 T=tf−tiT = t_f - t_i 均匀分割成 NN 个微小的时间步 Δt\Delta t

Δt=tf−tiN\Delta t = \frac{t_f - t_i}{N}

令 tk=ti+kΔtt_k = t_i + k \Delta t（其中 k=0,1,…,Nk = 0, 1, \dots, N），则 t0=tit_0 = t_i，tN=tft_N = t_f。于是利用指数算符的性质，时间演化算符可以拆解为 NN 个微小时间步演化算符的连乘

e−iℏH^(tf−ti)=(e−iℏH^Δt)N=e−iℏH^Δt⋅e−iℏH^Δt⋯e−iℏH^Δt⏟N 个e^{-\frac{i}{\hbar} \hat{H} (t_f - t_i)} = \left( e^{-\frac{i}{\hbar} \hat{H} \Delta t} \right)^N = \underbrace{ e^{-\frac{i}{\hbar} \hat{H} \Delta t} \cdot e^{-\frac{i}{\hbar} \hat{H} \Delta t} \cdots e^{-\frac{i}{\hbar} \hat{H} \Delta t} }_{N \text{ 个}}

于是一个整体的变换被我们切割为了连续变换

K(xf,tf;xi,ti)=∫ ⁣⋯∫dx1dx2…dxN−1∏k=1N⟨xk∣e−iℏH^Δt∣xk−1⟩K(x_f, t_f; x_i, t_i) = \int \dots \int dx_1 dx_2 \dots dx_{N-1} \prod_{k=1}^{N} \langle x_k \vert{} e^{-\frac{i}{\hbar} \hat{H} \Delta t} \vert{} x_{k-1} \rangle

这件事就可以等价地理解为对所有可能的（折线）路径求和，也就是费曼的路径积分思想。

下面计算某一特定的微小时间步

Uk,k−1=⟨xk∣e−iℏ(p^22m+V(x^))Δt∣xk−1⟩U_{k, k-1} = \langle x_k \vert{} e^{-\frac{i}{\hbar} \left( \frac{\hat{p}^2}{2m} + V(\hat{x}) \right) \Delta t} \vert{} x_{k-1} \rangle

使用Trotter 乘积公式并略去小量

e−iℏ(T^+V^)Δt=e−iℏT^Δte−iℏV^Δte^{-\frac{i}{\hbar}(\hat{T}+\hat{V})\Delta t} = e^{-\frac{i}{\hbar}\hat{T}\Delta t} e^{-\frac{i}{\hbar}\hat{V}\Delta t}
Uk,k−1≈⟨xk∣e−iℏp^22mΔte−iℏV(x^)Δt∣xk−1⟩U_{k, k-1} \approx \langle x_k \vert{} e^{-\frac{i}{\hbar}\frac{\hat{p}^2}{2m}\Delta t} e^{-\frac{i}{\hbar}V(\hat{x})\Delta t} \vert{} x_{k-1} \rangle

在其中插入动量本征态的完备性关系 ∫dpk∣pk⟩⟨pk∣=I^\int dp_k \vert{}p_k\rangle \langle p_k\vert{} = \hat{I}，并利用平面波内积 ⟨x∣p⟩=12πℏeiℏpx\langle x \vert{} p \rangle = \frac{1}{\sqrt{2\pi\hbar}} e^{\frac{i}{\hbar}px}

⟨xk∣e−iℏp^22mΔt∣xk−1⟩=∫dpk⟨xk∣pk⟩⟨pk∣e−iℏp^22mΔt∣xk−1⟩=∫dpk(12πℏeiℏpkxk)e−iℏpk22mΔt(12πℏe−iℏpkxk−1)=12πℏ∫dpkexp⁡[iℏ(pk(xk−xk−1)−pk22mΔt)]\begin{aligned} \langle x_k \vert{} e^{-\frac{i}{\hbar}\frac{\hat{p}^2}{2m}\Delta t} \vert{} x_{k-1} \rangle &= \int dp_k \langle x_k \vert{} p_k \rangle \langle p_k \vert{} e^{-\frac{i}{\hbar}\frac{\hat{p}^2}{2m}\Delta t} \vert{} x_{k-1} \rangle \\ &= \int dp_k \left( \frac{1}{\sqrt{2\pi\hbar}} e^{\frac{i}{\hbar}p_k x_k} \right) e^{-\frac{i}{\hbar}\frac{p_k^2}{2m}\Delta t} \left( \frac{1}{\sqrt{2\pi\hbar}} e^{-\frac{i}{\hbar}p_k x_{k-1}} \right) \\ &= \frac{1}{2\pi\hbar} \int dp_k \exp\left[ \frac{i}{\hbar} \left( p_k(x_k - x_{k-1}) - \frac{p_k^2}{2m}\Delta t \right) \right] \end{aligned}

求解动量积分后得到

⟨xk∣e−iℏp^22mΔt∣xk−1⟩=m2πiℏΔtexp⁡[iℏm(xk−xk−1)22Δt]\langle x_k \vert{} e^{-\frac{i}{\hbar}\frac{\hat{p}^2}{2m}\Delta t} \vert{} x_{k-1} \rangle = \sqrt{\frac{m}{2\pi i \hbar \Delta t}} \exp\left[ \frac{i}{\hbar} \frac{m (x_k - x_{k-1})^2}{2\Delta t} \right]

代回Uk,k−1U_{k, k-1}

Uk,k−1=m2πiℏΔtexp⁡{iℏ[m2(xk−xk−1Δt)2−V(xk−1)]Δt}U_{k, k-1} = \sqrt{\frac{m}{2\pi i \hbar \Delta t}} \exp\left\{ \frac{i}{\hbar} \left[ \frac{m}{2}\left(\frac{x_k - x_{k-1}}{\Delta t}\right)^2 - V(x_{k-1}) \right] \Delta t \right\}

会注意到指数的部分可以写成拉格朗日量的形式

将所有 NN 个时间步的积分连乘重新拼合回总传播子 K(xf,tf;xi,ti)K(x_f, t_f; x_i, t_i) 中

K(xf,tf;xi,ti)=lim⁡N→∞(m2πiℏΔt)N/2∫dx1…dxN−1exp⁡[iℏ∑k=1NL(xk,xk−xk−1Δt)Δt]K(x_f, t_f; x_i, t_i) = \lim_{N \to \infty} \left( \frac{m}{2\pi i \hbar \Delta t} \right)^{N/2} \int dx_1 \dots dx_{N-1} \exp\left[ \frac{i}{\hbar} \sum_{k=1}^N L\left(x_k, \frac{x_k - x_{k-1}}{\Delta t}\right) \Delta t \right]

求和在极限下变成积分，也即作用量

lim⁡Δt→0∑k=1NL(xk,x˙k)Δt=∫titfL(x(t),x˙(t)) dt=S[x(t)]\lim_{\Delta t \to 0} \sum_{k=1}^N L(x_k, \dot{x}_k) \Delta t = \int_{t_i}^{t_f} L(x(t), \dot{x}(t)) \, dt = S[x(t)]

把归一化常数与无限维的中间点积分统一写成路径测度 D[x(t)]\mathcal{D}[x(t)]

D[x(t)]≡lim⁡N→∞(m2πiℏΔt)N/2∏k=1N−1dxk\mathcal{D}[x(t)] \equiv \lim_{N \to \infty} \left( \frac{m}{2\pi i \hbar \Delta t} \right)^{N/2} \prod_{k=1}^{N-1} dx_k

最终得到费曼路径积分的核心表达式

K(xf,tf;xi,ti)=∫x(ti)=xix(tf)=xfD[x(t)] exp⁡(iℏS[x(t)])K(x_f, t_f; x_i, t_i) = \int_{x(t_i)=x_i}^{x(t_f)=x_f} \mathcal{D}[x(t)] \, \exp\left( \frac{i}{\hbar} S[x(t)] \right)

路径积分蒙特卡洛方法 / Path Integral Monte Carlo · 图 1

这个公式更多的或许是提供一个新的理解波函数传播的insight，毕竟对所有路径求和是一件很困难的事情。

2. PIMC

那么接下来你是不是觉得PIMC就是用蒙特卡洛来数值求解传播子呢？哼哼，很遗憾，并非如此。正如一开始说的，PIMC解决的问题与前面提到的波函数传播问题其实关系不大。

其真正的出发点是量子统计问题。我们知道，系统在某个温度 TT 下的配分函数 ZZ 以及任意物理量 A^\hat{A} 的热力学平均值 ⟨A⟩\langle A \rangle可以写作

Z=Tr(e−βH^)Z = \text{Tr}\left( e^{-\beta \hat{H}} \right)
⟨A⟩=1ZTr(e−βH^A^)\langle A \rangle = \frac{1}{Z} \text{Tr}\left( e^{-\beta \hat{H}} \hat{A} \right)

计算的难点其实就是e−βH^e^{-\beta \hat{H}}这一项，但聪明的你一下子发现这和刚才解决的问题如出一辙。等价地，将β\beta等价于虚时间 τ=ℏβ\tau = \hbar \beta，均匀分割为 MM 个片

τ=Δτ⋅M  ⟹  Δτ=βM\tau = \Delta \tau \cdot M \quad \implies \quad \Delta \tau = \frac{\beta}{M}

我们直接故技重施，把配分函数 ZZ 可以写成 MM 个因子连乘的迹

Z=Tr[(e−Δτ(T^+V^))M]Z = \text{Tr}\left[ \left( e^{-\Delta \tau (\hat{T} + \hat{V})} \right)^M \right]

求迹定义为

Tr(e−βH^)=∑n⟨n∣e−βH^∣n⟩\text{Tr}(e^{-\beta \hat{H}}) = \sum_{n} \langle n \vert{} e^{-\beta \hat{H}} \vert{} n \rangle

这里不需要是本征态，任意正交基底即可

于是这就类似于初末状态相同的传播子（相当于使用位置本征态），也即

Z=∫dx1∫dx2 ⁣⋯∫dxM ⟨x1∣e−ΔτH^∣xM⟩⟨xM∣e−ΔτH^∣xM−1⟩…⟨x2∣e−ΔτH^∣x1⟩Z = \int dx_1 \int dx_2 \dots \int dx_M \, \langle x_1 \vert{} e^{-\Delta \tau \hat{H}} \vert{} x_M \rangle \langle x_M \vert{} e^{-\Delta \tau \hat{H}} \vert{} x_{M-1} \rangle \dots \langle x_2 \vert{} e^{-\Delta \tau \hat{H}} \vert{} x_1 \rangle

且满足周期性边界条件：x0≡xMx_0 \equiv x_M

相似的计算过程就不写了，总之可以得到单步虚时间传播子

K(xk+1,xk;Δτ)≈(m2πℏ2Δτ)1/2exp⁡[−m2ℏ2Δτ(xk+1−xk)2−ΔτV(xk)]K(x_{k+1}, x_k; \Delta \tau) \approx \left( \frac{m}{2\pi \hbar^2 \Delta \tau} \right)^{1/2} \exp\left[ -\frac{m}{2\hbar^2 \Delta \tau} (x_{k+1} - x_k)^2 - \Delta \tau V(x_k) \right]

将 MM 个单步传播子拼合回配分函数 ZZ 中

Z=(m2πℏ2Δτ)M/2∫dx1dx2…dxMexp⁡(−∑k=1M[m2ℏ2Δτ(xk+1−xk)2+ΔτV(xk)])Z = \left( \frac{m}{2\pi \hbar^2 \Delta \tau} \right)^{M/2} \int dx_1 dx_2 \dots dx_M \exp\left( -\sum_{k=1}^M \left[ \frac{m}{2\hbar^2 \Delta \tau} (x_{k+1} - x_k)^2 + \Delta \tau V(x_k) \right] \right)

不妨写作

Z=C∫dx1…dxMe−βUeff(x1,…,xM)Z = C \int dx_1 \dots dx_M e^{-\beta U_{eff}(x_1, \dots, x_M)}

其中，有效势能函数定义为

Ueff(x1,…,xM)=∑k=1M[12kspring(xk+1−xk)2+1MV(xk)]U_{eff}(x_1, \dots, x_M) = \sum_{k=1}^M \left[ \frac{1}{2} k_{spring} (x_{k+1} - x_k)^2 + \frac{1}{M} V(x_k) \right]
kspring=mMℏ2β2k_{spring} = \frac{m M}{\hbar^2 \beta^2}

数学上可以视作是一个由 MM 个小球通过谐振子弹簧相连构成的环状聚合物

路径积分蒙特卡洛方法 / Path Integral Monte Carlo · 图 2

观察配分函数的形式，会发现概率恒正，我们可以试着进行蒙特卡洛采样。（相反，传播子问题由于指数上是个虚数，用蒙特卡洛就不太现实）

不妨直接以最简单的位置算符为示例，根据连续位置表象下求迹的直接公式

Tr(x^e−βH^)=∫−∞∞dx1 ⟨x1∣x^e−βH^∣x1⟩=∫−∞∞dx1 x1⋅⟨x1∣e−βH^∣x1⟩\text{Tr}(\hat{x} e^{-\beta \hat{H}}) = \int_{-\infty}^{\infty} dx_1 \, \langle x_1 \vert{} \hat{x} e^{-\beta \hat{H}} \vert{} x_1 \rangle = \int_{-\infty}^{\infty} dx_1 \, x_1 \cdot \langle x_1 \vert{} e^{-\beta \hat{H}} \vert{} x_1 \rangle

同样的用Trotter公式进行切分

Tr(x^e−βH^)=C∫dx1dx2…dxM (x1⋅e−βUeff(X))\text{Tr}(\hat{x} e^{-\beta \hat{H}}) = C \int dx_1 dx_2 \dots dx_M \, \Big( x_1 \cdot e^{-\beta U_{\text{eff}}(\mathbf{X})} \Big)

于是现在所求的物理量可以写作

⟨x⟩=Tr(x^e−βH^)Tr(e−βH^)=∫dX x1e−βUeff(X)∫dX e−βUeff(X)\langle x \rangle = \frac{\text{Tr}(\hat{x} e^{-\beta \hat{H}})}{\text{Tr}(e^{-\beta \hat{H}})} = \frac{\int d\mathbf{X} \, x_1 e^{-\beta U_{\text{eff}}(\mathbf{X})}}{\int d\mathbf{X} \, e^{-\beta U_{\text{eff}}(\mathbf{X})}}

采样时只有有限条路径，积分变求和

⟨x⟩=∑x1e−βUeff(X)∑e−βUeff(X)\langle x \rangle = \frac{\sum x_1 e^{-\beta U_{\text{eff}}(\mathbf{X})}}{\sum e^{-\beta U_{\text{eff}}(\mathbf{X})}}

注意到我们选取的是闭合路径，也就是说起点的选择满足置换对称性，不同的起点理应对应相同的权重，也就是说我们可以把起点的x1x_1换成更有对称性的xest(X)=1M∑k=1Mxkx_{\text{est}}(\mathbf{X}) = \frac{1}{M} \sum_{k=1}^{M} x_k（本质上是把一条路径换成了M条路径，因为具有相同权重所以在路径很多时肯定会收敛成同一个值，但在有限路径下具有统计优越性）

那么接下来做的事就是很典型的蒙特卡洛算法，定义路径构型 X\mathbf{X} 的概率密度函数为 P(X)=e−βUeff(X)ZP(\mathbf{X}) = \frac{e^{-\beta U_{\text{eff}}(\mathbf{X})}}{Z}。用 Metropolis 算法生成 NsampleN_{\text{sample}} 条符合 P(X)P(\mathbf{X}) 分布的环形路径样本

X(1),X(2),…,X(Nsample)\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \dots, \mathbf{X}^{(N_{\text{sample}})}

代入求和取平均，即为我们所求的统计意义上的位置量

⟨x⟩≈1Nsample∑a=1Nsample(1M∑k=1Mxk(a))\langle x \rangle \approx \frac{1}{N_{\text{sample}}} \sum_{a=1}^{N_{\text{sample}}} \left( \frac{1}{M} \sum_{k=1}^{M} x_k^{(a)} \right)

同理可以用来计算与位置直接相关的势能值

⟨V(x)⟩≈1Nsample∑a=1Nsample(1M∑k=1MV(xk)(a))\langle V(x) \rangle \approx \frac{1}{N_{\text{sample}}} \sum_{a=1}^{N_{\text{sample}}} \left( \frac{1}{M} \sum_{k=1}^{M} V(x_k)^{(a)} \right)

对非位置对角化的算子就比较麻烦，例如我们要求能量的统计平均值，可以利用公式

⟨E⟩=−∂ln⁡Z∂β=−1Z∂Z∂β\langle E \rangle = -\frac{\partial \ln Z}{\partial \beta} = -\frac{1}{Z} \frac{\partial Z}{\partial \beta}

写出离散化后的 Trotter-Suzuki 配分函数

Z=CM∫dXexp⁡[−S(X,β)]Z = C^M \int d\mathbf{X} \exp\left[ -S(\mathbf{X}, \beta) \right]
C=m2πℏ2Δτ=(mM2πℏ2β)1/2C = \sqrt{\frac{m}{2\pi \hbar^2 \Delta \tau}} = \left( \frac{m M}{2\pi \hbar^2 \beta} \right)^{1/2}
S(X,β)=∑k=1M[mM2ℏ2β(xk+1−xk)2+βMV(xk)]S(\mathbf{X}, \beta) = \sum_{k=1}^M \left[ \frac{m M}{2\hbar^2 \beta} (x_{k+1} - x_k)^2 + \frac{\beta}{M} V(x_k) \right]

带入求导可得

∂Z∂β=(∂CM∂β)∫dX e−S+CM∫dX (−∂S∂β)e−S\frac{\partial Z}{\partial \beta} = \left( \frac{\partial C^M}{\partial \beta} \right) \int d\mathbf{X} \, e^{-S} + C^M \int d\mathbf{X} \, \left( -\frac{\partial S}{\partial \beta} \right) e^{-S}
∂CM∂β=−M2βCM\frac{\partial C^M}{\partial \beta} = -\frac{M}{2\beta} C^M
∂S∂β=∂∂β∑k=1M[mM2ℏ2β(xk+1−xk)2+βMV(xk)]\frac{\partial S}{\partial \beta} = \frac{\partial}{\partial \beta} \sum_{k=1}^M \left[ \frac{m M}{2\hbar^2 \beta} (x_{k+1} - x_k)^2 + \frac{\beta}{M} V(x_k) \right]
=−∑k=1MmM2ℏ2β2(xk+1−xk)2+1M∑k=1MV(xk)= -\sum_{k=1}^M \frac{m M}{2\hbar^2 \beta^2} (x_{k+1} - x_k)^2 + \frac{1}{M} \sum_{k=1}^M V(x_k)
⟨E⟩=∫dX P(X)[M2β−∑k=1MmM2ℏ2β2(xk+1−xk)2+1M∑k=1MV(xk)]⏟Ethermo(X)\langle E \rangle = \int d\mathbf{X} \, P(\mathbf{X}) \underbrace{\left[ \frac{M}{2\beta} - \sum_{k=1}^M \frac{m M}{2\hbar^2 \beta^2} (x_{k+1} - x_k)^2 + \frac{1}{M} \sum_{k=1}^M V(x_k) \right]}_{E_{\text{thermo}}(\mathbf{X})}

于是即有能量对应的可用于蒙特卡洛计算的估计算子Ethermo(X)E_{\text{thermo}}(\mathbf{X})，其余的一些热力学相关量也可用类似方法导出。

3. 多体推广与负符号问题

以上讨论的是单粒子的系综平均问题，但这一思路很容易地可以扩展到多粒子体系。

推广到 NN 个粒子的多体系统时，只需要 N×MN\times M 的矩阵 R=(r1,…,rN)\mathbf{R} = (\mathbf{r}_1, \dots, \mathbf{r}_N)，其中每个时间切片上包含所有 NN 个粒子的位置（为了方便还是接着在一维上讨论）

Ueff(R)=∑k=1M[∑i=1NmM2ℏ2β2(ri,k+1−ri,k)2+1MV(r1,k,…,rN,k)]U_{\text{eff}}(\mathbf{R}) = \sum_{k=1}^M \left[ \sum_{i=1}^N \frac{m M}{2\hbar^2 \beta^2} (\mathbf{r}_{i, k+1} - \mathbf{r}_{i, k})^2 + \frac{1}{M} V(\mathbf{r}_{1, k}, \dots, \mathbf{r}_{N, k}) \right]

对于 NN 个粒子在逆温度 β=1/kBT\beta = 1/k_B T 下的量子多体系统，其正则配分函数本质上依然是对密度算符求迹

Z=Tr(e−βH^)Z = \text{Tr}\left( e^{-\beta \hat{H}} \right)

但是多体系统的交换对称性和反对称性使得其连续位置基底有一些差异，具体一点说，假设我们先不考虑粒子全同性，则完备的位置基底是单粒子空间的直积

∣r⟩≡∣r1⟩(1)⊗∣r2⟩(2)⊗⋯⊗∣rN⟩(N)\vert{}\mathbf{r}\rangle \equiv \vert{}\mathbf{r}_1\rangle^{(1)} \otimes \vert{}\mathbf{r}_2\rangle^{(2)} \otimes \dots \otimes \vert{}\mathbf{r}_N\rangle^{(N)}

满足正交性 ⟨r∣r′⟩=∏i=1Nδ(ri−ri′)=δN(r−r′)\langle \mathbf{r} \vert{} \mathbf{r}' \rangle = \prod_{i=1}^N \delta(\mathbf{r}_i - \mathbf{r}_i') = \delta^N(\mathbf{r} - \mathbf{r}') 和完备性 ∫dr1…drN ∣r⟩⟨r∣=I^dist\int d\mathbf{r}_1 \dots d\mathbf{r}_N \, \vert{}\mathbf{r}\rangle \langle \mathbf{r}\vert{} = \hat{I}_{\text{dist}}

真实体系的全同性要求系统的物理状态必须处于置换群 SNS_N 的一维表示，定义置换算符 P^σ\hat{P}_\sigma 对坐标基底的作用（σ∈SN\sigma \in S_N）

P^σ∣r1,…,rN⟩=∣rσ−1(1),…,rσ−1(N)⟩≡∣Pσr⟩\hat{P}_\sigma \vert{}\mathbf{r}_1, \dots, \mathbf{r}_N\rangle = \vert{}\mathbf{r}_{\sigma^{-1}(1)}, \dots, \mathbf{r}_{\sigma^{-1}(N)}\rangle \equiv \vert{}P_\sigma \mathbf{r}\rangle

由此，我们可以定义全同粒子投影算符 P^±\hat{\mathcal{P}}_\pm，正负号分别代表交互对称性和反对称性

P^±=1N!∑σ∈SN(±1)sign(σ)P^σ\hat{\mathcal{P}}_\pm = \frac{1}{N!} \sum_{\sigma \in S_N} (\pm 1)^{\text{sign}(\sigma)} \hat{P}_\sigma

将投影算符作用于未对称化基底 ∣r⟩\vert{}\mathbf{r}\rangle，可以构造出全同粒子的连续位置态 ∣r⟩±\vert{}\mathbf{r}\rangle_\pm

∣r⟩±≡N! P^±∣r⟩=1N!∑σ∈SN(±1)sign(σ)∣Pσr⟩\vert{}\mathbf{r}\rangle_\pm \equiv \sqrt{N!} \, \hat{\mathcal{P}}_\pm \vert{}\mathbf{r}\rangle = \frac{1}{\sqrt{N!}} \sum_{\sigma \in S_N} (\pm 1)^{\text{sign}(\sigma)} \vert{}P_\sigma \mathbf{r}\rangle

会发现满足置换对称性的位置态都被我们视作了同一个态。

这组基底满足比较特殊的正交完备性，我们可以试着计算两个全同粒子位置态 ∣r⟩±\vert{}\mathbf{r}\rangle_\pm 与 ∣r′⟩±\vert{}\mathbf{r}'\rangle_\pm 的内积

±⟨r∣r′⟩±=1N!∑σ∈SN∑τ∈SN(±1)sign(σ)+sign(τ)⟨Pσr∣Pτr′⟩=1N!∑σ∈SN∑τ∈SN(±1)sign(σ−1τ)⟨r∣Pσ−1τr′⟩\begin{aligned} {}_\pm\langle \mathbf{r} \vert{} \mathbf{r}' \rangle_\pm &= \frac{1}{N!} \sum_{\sigma \in S_N} \sum_{\tau \in S_N} (\pm 1)^{\text{sign}(\sigma) + \text{sign}(\tau)} \langle P_\sigma \mathbf{r} \vert{} P_\tau \mathbf{r}' \rangle \\ &= \frac{1}{N!} \sum_{\sigma \in S_N} \sum_{\tau \in S_N} (\pm 1)^{\text{sign}(\sigma^{-1}\tau)} \langle \mathbf{r} \vert{} P_{\sigma^{-1}\tau} \mathbf{r}' \rangle \end{aligned}

令 μ=σ−1τ\mu = \sigma^{-1}\tau，当 σ\sigma 固定时，对 τ\tau 的求和等价于对 μ∈SN\mu \in S_N 的求和

±⟨r∣r′⟩±=∑μ∈SN(±1)sign(μ)δN(r−Pμr′){}_\pm\langle \mathbf{r} \vert{} \mathbf{r}' \rangle_\pm = \sum_{\mu \in S_N} (\pm 1)^{\text{sign}(\mu)} \delta^N(\mathbf{r} - P_\mu \mathbf{r}')

这其实暗示我们全同粒子的位置内积不仅包括 r=r′\mathbf{r} = \mathbf{r}' 的直通项，还包括将 r′\mathbf{r}' 中的粒子位置进行任意置换后与 r\mathbf{r} 重合的微观交换项。

至于完备性，我们将全同粒子连续位置基底 ∣r⟩±\vert{}\mathbf{r}\rangle_\pm 对全空间积分

∫dr ∣r⟩±±⟨r∣=∫dr (N!P^±∣r⟩)(N!⟨r∣P^±)=N! P^±(∫dr ∣r⟩⟨r∣)P^±=N! P^±I^distP^±=N! I^±\begin{aligned} \int d\mathbf{r} \, \vert{}\mathbf{r}\rangle_\pm {}_\pm\langle \mathbf{r}\vert{} &= \int d\mathbf{r} \, \left( \sqrt{N!} \hat{\mathcal{P}}_\pm \vert{}\mathbf{r}\rangle \right) \left( \sqrt{N!} \langle \mathbf{r}\vert{} \hat{\mathcal{P}}_\pm \right) \\ &= N! \, \hat{\mathcal{P}}_\pm \left( \int d\mathbf{r} \, \vert{}\mathbf{r}\rangle \langle \mathbf{r}\vert{} \right) \hat{\mathcal{P}}_\pm \\ &= N! \, \hat{\mathcal{P}}_\pm \hat{I}_{\text{dist}} \hat{\mathcal{P}}_\pm = N! \, \hat{I}_\pm \end{aligned}

两边同时除以 N!N!，即得到全同粒子连续位置完备性关系

1N!∫dr1…drN ∣r⟩±±⟨r∣=I^±\frac{1}{N!} \int d\mathbf{r}_1 \dots d\mathbf{r}_N \, \vert{}\mathbf{r}\rangle_\pm {}_\pm\langle \mathbf{r}\vert{} = \hat{I}_\pm

有了这些理论准备后，我们可以开始推导多体系统的配分函数形式了，同样从最基础的求迹公式出发，利用位置基底进行展开

Z=TrH±(e−βH^)Z = \text{Tr}_{\mathcal{H}_\pm} \left( e^{-\beta \hat{H}} \right)
Z±=1N!∫dr ±⟨r∣e−βH^∣r⟩±Z_\pm = \frac{1}{N!} \int d\mathbf{r} \, {}_\pm\langle \mathbf{r}\vert{} e^{-\beta \hat{H}} \vert{}\mathbf{r}\rangle_\pm

然后就是Trotter切分和路径积分这些比较繁琐的过程，具体推导就略去吧，我们直接看结果

对于玻色子，取 (±1)sign(σ)=+1(\pm 1)^{\text{sign}(\sigma)} = +1

ZBoson=1N!∑σ∈SN∫dR1…dRM exp⁡[−Seff(R;σ)]Z_{\text{Boson}} = \frac{1}{N!} \sum_{\sigma \in S_N} \int d\mathbf{R}_1 \dots d\mathbf{R}_M \, \exp\left[ -S_{\text{eff}}(\mathbf{R}; \sigma) \right]

对费米子，取 (±1)sign(σ)=(−1)sign(σ)(\pm 1)^{\text{sign}(\sigma)} = (-1)^{\text{sign}(\sigma)}

ZFermion=1N!∑σ∈SN(−1)sign(σ)∫dR1…dRM exp⁡[−Seff(R;σ)]Z_{\text{Fermion}} = \frac{1}{N!} \sum_{\sigma \in S_N} \mathbf{(-1)^{\text{sign}(\sigma)}} \int d\mathbf{R}_1 \dots d\mathbf{R}_M \, \exp\left[ -S_{\text{eff}}(\mathbf{R}; \sigma) \right]

其中

Seff(R;σ)=∑k=1M[∑i=1NmM2ℏ2β2(ri,k+1−ri,k)2+βMV(r1,k,…,rN,k)]S_{\text{eff}}(\mathbf{R}; \sigma) = \sum_{k=1}^M \left[ \sum_{i=1}^N \frac{m M}{2\hbar^2 \beta^2} (\mathbf{r}_{i, k+1} - \mathbf{r}_{i, k})^2 + \frac{\beta}{M} V(\mathbf{r}_{1, k}, \dots, \mathbf{r}_{N, k}) \right]

其连续边界与原先的直接取等不同，可以由置换 σ\sigma 确定

ri,M+1=rσ(i),1\mathbf{r}_{i, M+1} = \mathbf{r}_{\sigma(i), 1}

然后就是蒙特卡洛的采样过程，但这里有一个令人在意的点，对于费米子，原先在单粒子体系中某一路径构型对应的恒正的概率权重W(C)W(\mathcal{C})在此处变成了

W(C)Fermion=(−1)sign(σ)e−Seff(R;σ)W(\mathcal{C})_{\text{Fermion}} = (-1)^{\text{sign}(\sigma)} e^{-S_{\text{eff}}(\mathbf{R}; \sigma)}

这个东西是正负不定的，与置换的奇偶有关。这意味着我们没法将其等价地视作概率分布，毕竟根据路径积分，我们需要计算的物理量是

⟨A⟩Fermion=∫dC A(C)⋅WFermion(C)∫dC WFermion(C)\langle A \rangle_{\text{Fermion}} = \frac{\int d\mathcal{C} \, A(\mathcal{C}) \cdot W_{\text{Fermion}}(\mathcal{C})}{\int d\mathcal{C} \, W_{\text{Fermion}}(\mathcal{C})}

不过形式上，我们仍然可以将任意构型的权重 WFermion(C)W_{\text{Fermion}}(\mathcal{C}) 拆分为符号与绝对值的乘积

WFermion(C)=sign(C)⋅∣WFermion(C)∣W_{\text{Fermion}}(\mathcal{C}) = \text{sign}(\mathcal{C}) \cdot \vert{}W_{\text{Fermion}}(\mathcal{C})\vert{}

定义一个由绝对值权重构成的参考辅助配分函数 ZrefZ_{\text{ref}}

Zref=∫dC ∣WFermion(C)∣>0Z_{\text{ref}} = \int d\mathcal{C} \, \vert{}W_{\text{Fermion}}(\mathcal{C})\vert{} > 0

于是得到一个形式上等价的方程

⟨A⟩Fermion=1Zref∫dC A(C)⋅sign(C)⋅∣WFermion(C)∣1Zref∫dC sign(C)⋅∣WFermion(C)∣\langle A \rangle_{\text{Fermion}} = \frac{\frac{1}{Z_{\text{ref}}} \int d\mathcal{C} \, A(\mathcal{C}) \cdot \text{sign}(\mathcal{C}) \cdot \vert{}W_{\text{Fermion}}(\mathcal{C})\vert{}}{\frac{1}{Z_{\text{ref}}} \int d\mathcal{C} \, \text{sign}(\mathcal{C}) \cdot \vert{}W_{\text{Fermion}}(\mathcal{C})\vert{}}

因为 ∣WFermion(C)∣>0\vert{}W_{\text{Fermion}}(\mathcal{C})\vert{} > 0，所以 Pref(C)≡∣WFermion(C)∣ZrefP_{\text{ref}}(\mathcal{C}) \equiv \frac{\vert{}W_{\text{Fermion}}(\mathcal{C})\vert{}}{Z_{\text{ref}}} 完全满足概率测度的所有条件，我们依然可以使用蒙特卡洛方法进行采样。

换句话说我们强行转移符号，把概率的部分提取出来。但这么做的代价就是会产生严重的正负抵消问题，采样数必须指数增长才能确保结果的精度，也即大名鼎鼎的负符号问题，与DMC遇到的问题如出一辙。

关于PIMC最出名的成果应该是Ceperley关于液态4He^4\text{He}的超流相变相关的工作Path integrals in the theory of condensed helium，希望自己有空时能去认真研读一下。

那么，以上。

继续阅读 · 计算物理

HF Method ＆ Post HF
← 返回笔记目录回到顶部 ↑