Yixiong's Blog
0%

Feynman path integrals

Posted on Edited on In
Symbols count in article: 2k Reading time ≈ 4 mins.

考虑到状态的演化，最直接的想法是通过经典力学通过拉格朗日量来构建经典的运动路径。如何处理的问题不是精确的状态，而是概率密度的分布，如何求概率密度分布的演化？路径积分则是十分自然的想法。将每一个小的状态，以任意的路径进行演化，最后得到最终转态的概率分布。这样处理的麻烦在于，不是每一条路径都是等价的，有一些路径发生的概率高一些，另一些则低一些。如何表示这样的概率分布，以及实际得到演化后的概率分布，则是本文要探讨的。

当然，这只是一个简单的版本，需要更多的讨论。

参考：

路径积分的基本构架

本部分的目标： 1. 构建传播子 2. 相空间、位形空间的形式 3. 以自由粒子为例，展示两种形式

传播子是研究态随时间演化的方法，这样演化过程由基本方程确定，在量子力学中这样的方程是薛定谔方程。但是随时间演化的方程有很多，这些同样可以在传播子的框架下研究。

由于薛定谔方程过于出名，接下来的讨论在薛定谔方程下进行。

传播子

定义一个算符 $U$，它的作用是给出演化结果:

然后求解 $U$ 的表达式，研究时间演化最熟悉的是薛定谔方程:

如果时间不变自然不发生演化，以此作为初始条件，即 $t=t_0 \Rightarrow U=I$。确定系数之后得到了么正时间演化算符的表达式:

在形式上地使用：

接下来$\eqref{1_U}$式取坐标表象，即左乘 $\langle\vec{r}|$，再插入完备性关系式 $I=\int\left|\vec{r}_0\right\rangle\left\langle\vec{r}_0\right| \mathrm{d} \vec{r}_0$，可得：

其中 $K\left(\vec{r}, t ; \vec{r}_0, t_0\right)$ 称作传播子, 作用如上所示, 路径积分的最终目的就是得到其具体表达式。

传播子的性质

设置时间节点（无特殊声明全文适用）：

注意不平均分割时间。由传播子的定义，我们不难发现它具有如下特性：

注意我们要求上式中恒满足关系 $\vec{r}_N=\vec{r}, t_N=t, a>b \Leftrightarrow t_a>t_b$。利用该性质进行无穷分割:

将积分号里边儿的那个无穷小过程传播子换一个符号, 记为:

相空间的路径积分

要求 $K\left(\vec{r}, t ; \vec{r}_0, t_0\right)$ ，需先处理 $\mathcal{K}\left(\vec{r}_n, t_n ; \vec{r}_{n-1}, t_{n-1}\right)$ 。

$\mathcal{K}$表示极短时间内从一个状态到另一个状态，在极限情况下式精确的；$K$是一段时间的跳跃。因此有 $K \neq \mathcal{K}$，在一些情况下这两个是相等的例如自由粒子。

可以推知：

将上式代回传播子的表达式就得到相空间的路径积分:

其中的 $\mathcal{D}$ 表征对各种怪异路径进行积分。注意位置积分比动量积分少一重，因为端点是固定的。

位形空间的路径积分

位形空间的路径积分，就是把动量先积掉：

显然后面的积分在自由粒子的格林函数中求解过, 这里直接而代入结论:

将上式代回无穷小间隔传播子表达式得:

将上述结果代回传播子的表达式得到位形空间的路径积分:

传播子的两类表达式

综上所述, 路径积分中的传播子就表达为:

实际运算中需要用到:

其中 $\left\{\begin{array}{l}\Delta t_n=t_n-t_{n-1} \\ \Delta \vec{r}_n=\vec{r}_n-\vec{r}_{n-1} \\ \dot{\vec{r}}_n=\Delta \vec{r}_n / \Delta t_n \\ n \in\{1,2,3, \cdots, N\}\end{array} \quad\right.$ 且规定 $\left\{\begin{array}{l}\vec{r}_N=\vec{r} \\ t_N=t \\ t_0b \Leftrightarrow t_a>t_b \\ \min \left\{t_n-t_{n-1}\right\} \rightarrow 0\end{array}\right.$
根号下的虚数单位 $i$ 定义为 $\exp [i \pi / 2]$。

自由粒子 - 利用对角化求解路径积分

已完成内容： 1. 路径积分的两种形式 2. 自由粒子的路径积分 本部分的目标： 1. 讨论自由粒子在位形空间的路径积分，得到传播子的具体形式 2. 讨论自由粒子在相空间的路径积分，得到传播子的具体形式所属

对时间采用平均分割的方案$\Delta t_n=\frac{t-t_0}{N}=\varepsilon$，在极限情况下讨论$N\to \infty\quad \text{or}\quad \varepsilon\to\infty$。由$\eqref{1_weixin}$可知自由粒子位形空间传播子为：

将位置$x_n$表示成经典路径$x^c_n$与路径扰动$y_n$两个部分，其中下标表示时刻$t_n$：

并且有边界条件，在端点是固定不动的：

以及关系：

代入$\eqref{free_partial_S}$后：

分为三个部分：

Path Integral

如上图所示，表示在固定初始与终点的路径积分，每条线表示可能的路径，其中红线为经典路径。

针对交叉项：

考虑到自由粒子$F=ma=m\ddot{x}^{\mathrm{c}}=0$，因此交叉项$S^{\times}=0$。

对于经典路径项$S^{\mathrm{c}}$：

经典自由粒子解的 $\dot{x}^c$ 是一个常数, 且知道 $\dot{x}^{\mathrm{c}}=\frac{x-x_0}{t-t_0}$。所以：

端点给定，并且把经典路径分离出来，因此可以将传播子写为：

接下来要化简的是$S^{\prime}$，考虑到在两端点附近$y_0=y_N=0$：

其中$A$是一个实对称矩阵，$A_{mn}=2\delta_{m,n}-\delta_{m,n-1}-\delta_{m-1,n}$，如果不去掉端点的话，就会造成不是一个实对称矩阵。既然已经写成对称矩阵的形式了，那接下来就是精确对角化，将其耦合项解耦，从而方便积分。

其中 $R$ 是正交矩阵，满足 $R^{\mathrm{T}} R=|R|=1$ 同时注意么正变换不改变本征值:

令:

再进行变量代换:

因此传播子可以写为：

从$\eqref{3_propagator}$可以看出，传播子最后只与特征值的连乘积有关，下面就是求解$\prod_{n=1}^{N-1} \lambda_n$的值。求解本征的可通过本征多项式直接得到。根据$A$ 的本征多项式 $|A-\lambda I|=P_{N-1}(\lambda)$，构造多项式：

将$P_n(\lambda)$从第一行展开得到：

定义一个能提升多项式阶数的线性算子$L:LP_{n-1}=P_n$，然后将其代入递推关系$\eqref{2_P}$：

接下来分别求解$(L-\alpha)P_n=0$与$(L-\beta)P_n=0$：

同理可得：

因此通解为：

根据多项式$P_1=2-\lambda$与$P_2=(2-\lambda)^2-1$可知$P_1=1$，由此解得：

对于连乘积有$\prod_{n=1}^{N-1} \lambda_n=|A|$，根据$|A-\lambda I|=P_{N-1}(\lambda)$可得到$|A|=P_{N-1}(0)=N$，并将结论代入$\eqref{3_propagator}$：

谐振子 - 利用对角化求解路径积分

已完成内容： 1. 自由粒子的路径积分 本部分的目标： 1. 讨论谐振子的路径积分

对于一维谐振子势$V=\frac{1}{2}m\omega^2X^2$，研究从$t_0\to t$的演化结果。位形空间传播子：

同样进行代换，分离为经典路径和扰动项$x_n=x_n^{\mathrm{c}}+y_n$：

同样分为三部分：

先利用$y_n$的边界条件计算交叉项$S^\times$：

根据经典方程$m\ddot{x}^{\mathrm{c}}=-\frac{\partial V}{\partial x^{\mathrm{c}}}=-m\omega^2x^\mathrm{c}$，则有$S^\times=0$。

再计算经典路径$S^\mathrm{c}$：

接下来计算经典路径的形式：

1
2
DSolve[{m D[x[t], {t, 2}] == -m \[Omega]^2 x[t], x[t0] == x0, 
   x[tend] == xend}, x[t], t] // Simplify

得到结果：

将经典路径的表达形式$\eqref{classical}$代入经典路径的路径积分中$\eqref{3_classical}$，可以得到路径积分的具体表达式。

1
2
3
4
5
x[t] = -Csc[(t0 - tend) \[Omega]] (xend Sin[(t - t0) \[Omega]] - 
     x0 Sin[(t - tend) \[Omega]]);
Integrate[
  m/2 (D[x[t], {t, 1}]^2 - \[Omega]^2 x[t]^2), {t, t0, 
   tend}] // FullSimplify

路径 $x_n^{\mathrm{c}}=x^{\mathrm{c}}\left(t_n\right)$ 都是给定的,故:

接下来需要解决的是：

首先将$\Delta y_n$利用差分的办法改写为矩阵形式：

与自由粒子一样，因为是实正定矩阵，可以将其对角化。

接下来解决连乘积$\prod_{n=1}^{N-1}\left(\lambda_n-\varepsilon^2 \omega^2\right)$，对$\eqref{3_Sprime}$的本征多项式入手，与自由粒子的区别在于对角线上会多出来一项：

对第一行进行展开，可以得到递推关系：

定义线性算子$LP_{n}(\lambda)=P_{n+1}(\lambda)$，并将其代入递推关系可以得到：

进行分解：

接下来分别求解$(L-\beta)P_{n}=0$与$(L-\alpha)P_{n}=0$:

因此$P$通项写为：

代入初始条件：

由递推关系$\eqref{3_ditui}$可以得到$P_0=1$，$P_n$的表达式：

通过以上讨论已经成功将特征多项$\eqref{3_ditui}$通项表达形式写出来。接下来需要具体写出 $\alpha$ 与 $\beta$ 的表达形式。结合需要具体求解的问题$\prod_{n=1}^{N-1}\left(\lambda_n-\varepsilon^2 \omega^2\right)=|P_{N-1}(\lambda=0)|$，可知$\lambda=0$:

可以得到：

所以谐振子的传播子写为：

应该考虑相位了，但现在还没出现问题，之再说吧。详细讨论看上面引用的内容。

路径积分蒙特卡洛

参考：

考虑量子力学中的哈密顿量以及波函数：

在使用蒙卡数值计算的过程中存在两个问题：

  • 指数是复数，不能用使用正实数概率分布表示
  • $H$是量子力学算符，动能与势能不对易

Continuing to Imaginary Time

解决复数问题。

通过引入虚时$\tau=-it$，将复数变为一个实数处理，这种做法也出现在扩散蒙特卡洛（diffusion Monte Carlo, DMC）。通过定义配分函数，产生与统计力学的联系：

Discretizing the Time Dimension

解决不对易的问题。

将时间$t$放在$M$的时间格点上，$\Delta \tau = \tau/M$，这样有$M-1$个中间时间。在每一个中间的时间步骤中都有完整的本征态：

然后沿着时间的顺序将配分函数分解：

根据时间分割，将一段时间的演化过程，转变为几个“瞬间”的演化过程$\Delta \tau\to 0$。这其实就是路径积分的思想。

接下来通过这样的小量，解决不对易的关系。根据 Baker-Campbell-Hausdorff (BCH)定理：

把这个关系用在短时演化的过程中：

通过BCH定理可以知道，将$H$分为两项，需要计算交叉项。但是这个交叉项是$\Delta\tau$的高阶小量可以直接忽略。

Path Integral Monte Carlo Algorithm

配分函数的最终表达式为：

看这个形式，对比统计晶格中的形式。具有一定的对应关系，路径积分等价于晶格所有构型的求和，因此真正需要关心的是“能量项”：

将每一格运动区间进行计算，从而得出最终的“能量”，更准确应该是作用量。然后利用 Metropolis algorithm 进行重要性抽样，将不同权重的路径进行加权运算。

Code

以下是一个对于一维谐振子的实现，核心部分在sample模块中。这个实现，每次只对前一次路径中一步进行修改，然后接受概率设为exp(-DeltaE * Delta_t / hbar)，这里作者将$\hbar=1$进行无量纲化了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
class PathIntegral:
    def __init__(self, nt=100, tmax = 2*np.pi):
        self.nt = nt
        self.tmax = tmax
        self.t = np.linspace(0, self.tmax, nt)
        self.sample_idx = []
        
    def sample(self, steps=100000):
        # 1. initialize the path
        x = np.zeros_like(self.t)
        dt = self.t[1] - self.t[0]
        ilist = np.random.choice(range(1, self.nt-1), size=steps)
        dxlist = np.random.uniform(-1, 1, size=steps)
        sampled = []           
        # randomly pick x[i]
        # do random walk x[i] -> x[i] + r 
        # if DeltaE < 0: accept
        # elif DeltaE > 0: accept with probability exp(-DeltaE * Delta_t / hbar)
        for n in tqdm(range(steps)):
            i = ilist[n]
            dx = dxlist[n]
            
            #if x[i] + dx <= 0.0: continue
            ds = (self.kinetic_energy(x[i-1], x[i]+dx, dt) -
                  self.kinetic_energy(x[i-1], x[i], dt) +
                  self.kinetic_energy(x[i]+dx, x[i+1], dt) -
                  self.kinetic_energy(x[i], x[i+1], dt) + 
                  self.potential(x[i]+dx) -
                  self.potential(x[i])
                 ) * dt
            
            if ds < 0: 
                x[i] = x[i] + dx
            else:
                r = np.random.rand()
                if r < np.exp(-ds):
                    x[i] = x[i] + dx
            sampled.append(x[i])
            
            self.sample_idx.append(i)
        return sampled, self.sample_idx
    
    def kinetic_energy(self, x1, x2, dt):
        return 0.5 * (x2 - x1)**2 / dt**2
    
    def potential(self, x):
        return 0.5 * x**2
    
    def plot(self):
        tmax = self.tmax
        plt.vlines(self.t, ymin=-tmax, ymax=tmax, color='k', linestyles='dotted', alpha=0.9)
        dt = self.t[1] - self.t[0]
        xi = np.random.uniform(-1, 1, self.nt)
        xi[0] = 0
        xi = np.cumsum(xi)
        xi[-1] = 0
        plt.plot(self.t, xi)
        plt.ylabel(r'x $\rightarrow$', loc='center')
        plt.xlabel(r't $\rightarrow$', loc='center')