经典网文分享: 软物质的统计力学（三）

五、Fokker-Planck方程

之前的朗之万方程描述的是某个流体中悬浮着的粒子的轨迹，但一个随机过程的"导数"实在是一个麻烦重重含混不清的概念。现在我们终于可以回到一个让人舒适一点的框架中了。对于一个随机变量，它的概率密度函数已经包含了这个随机变量的全部信息。对于随机过程也是如此，只是这个分布函数是含时的了。我们记单粒子概率密度函数 $P(x;t)$ 是发现一t时刻随机变量取值在x附近的概率密度（分号是用来把概率密度的变量和参数分开。后文有时候会不严格遵守这种写法。）。或者使用条件概率，将初始条件也计入其中， $P(x;t|x_0;t_0)$ 为变量在 $t_0$ 时等于 $x_0$ 的条件下，t时刻取值在x附近的概率密度。那么只要能够计算这个概率密度函数，我们就知道了系统的一切信息。值得注意的是，我们悄悄换掉了变量x的涵义。在之前的讨论中它指某一个随机变量的具体取值，现在它指实轴上的某个位置。

我们现在依然在考虑无时滞的系统，所以将讨论依然限定在马尔可夫过程中。马尔可夫过程意味着如果考虑随机变量的一条轨迹 $(t_0,x_0),\ (t_1,x_1),\ \cdots,\ (t_n,x_n)$ ，那么联合概率密度 $P(x_n,t_n,x_{n-1},t_{n-1},\cdots,x_0,t_0)=P(x_n,t_n|x_{n-1},t_{n-1},\cdots,x_0,t_0)P(x_{n-1},t_{n-1}|x_{n-2},t_{n-2},\cdots,x_0,t_0)\cdots P(x_0,t_0)$

中的条件概率只与最近的时刻有关，即

$P(x_n,t_n,x_{n-1},t_{n-1},\cdots,x_0,t_0)=P(x_n,t_n|x_{n-1},t_{n-1})P(x_{n-1},t_{n-1}|x_{n-2},t_{n-2})\cdots P(x_0,t_0)$

即我们称为历史无关。于是对于任选的一个中间时刻 $t>t'>t_0$ ，穷尽随机变量在这时刻可能的取值，我们有Chapman-Komogorov方程

$P(x;t|x_0;t_0)=\int_\mathbb{R} P(x;t|x',t',x_0,t_0)dx'=\int_\mathbb{R} P(x;t|x';t')P(x';t'|x_0,t_0)dx'$

我们定义态转移函数 $W(x\to x')dt=P(x';t+dt|x;t)$ （如果 $x\neq x'$ ），为dt时间内变量从x附近变至x'附近的概率密度。显然dt时间后变量还留在x附近的概率密度是 $P(x;t+dt|x;t)=1-dt\int_{x\neq x'} W(x\to x')dx'$ 。套用CK方程即

$P(x;t+dt)=dt\int_{x\neq x'} W(x'\to x)P(x';t)dx'+P(x;t)\left(1-dt\int_{x\neq x'} W(x\to x')dx'\right)$

即（因为对称性，积分有没有去掉 $x=x'$ 反而不重要了）

$\frac{\partial P(x;t)}{\partial t}=\int W(x'\to x)P(x';t)dx'-\int W(x\to x')P(x;t)dx'$

这表示马尔可夫过程的概率密度随时间的增量等于这段时间内流入这一区域的概率，减去流出这一区域的概率（平时我们就是这样直接写出主方程的）。或者对于离散随机过程，

$\frac{\partial P(x;t)}{\partial t}=\sum_{x'}W(x'\to x)P(x';t)-\sum_{x'} W(x\to x')P(x;t)$

这个方程称为（正向）主方程，描述了概率密度函数随时间的变化。或者利用马尔可夫过程的时间平移对称性，类似地还有

$-\frac{\partial P(x;t|x_0;t_0)}{\partial t_0}=\int W(x_0\to x')(P(x;t|x';t_0)-P(x;t|x_0;t_0))dx'$

称为反向主方程。

重要的基本概念先介绍到这里，实际上，我们将要得到的Fokker-Planck方程，就是朗之万方程对应的马尔可夫过程的主方程。（我曾以为所有连续状态空间的主方程都是FP方程，但老板告诉我FP方程往往只指与朗之万方程对应的主方程，其他场合我们依然沿用主方程的称呼。）这里我们使用一个物理上比较直观但数学上不够严格的推导。假设一个粒子的轨迹 $x(t)$ （作为一个随机过程）满足朗之万方程（记号和以前维持一致，头上加点一般表示对时间求导）

$\dot{x}(t)=f(x,t)+\sqrt{2D(x,t)}\xi(t)$

我们考虑任意一个 $C^2$ 光滑的实函数 $g(x)$ ，我们定义它在t时刻的（系综）期望是

$\langle g(x(t)|x(0)=x_0)\rangle\equiv\int_\mathbb{R} dx'\, g(x')P(x';t|x_0;0)$

对时间求导，有

$\int_\mathbb{R} dx'\,g(x')\frac{\partial P(x';t|x_0,0)}{\partial t}=\left\langle\frac{d g(x(t)|x(0)=x_0)}{dt}\right\rangle$

这时候我们要用伊藤积分公式了，等式右边于是等于

$\left\langle\frac{d g(x(t)|x(0)=x_0)}{dt}\right\rangle=\left\langle\frac{d g}{dx}\frac{dx}{dt}+D(x,t)\frac{d^2g}{dx^2}\right\rangle=\left\langle\frac{d g}{dx}(f(x,t)+\sqrt{2D(x,t)}\xi(t))+D(x,t)\frac{d^2g}{dx^2}\right\rangle$

还记得伊藤积分的特性告诉我们 $\xi$ 和 $x$ 无关，而且因为 $\xi$ 是高斯分布，其系综平均为 $\langle\xi\rangle=0$ （你如果感觉到这里不太严格，很正常，因为 $\xi$ 本身不良定）。那么按照定义式，

$\left\langle\frac{dg}{dx}f(x,t)+D(x,t)\frac{d^2g}{dx^2}\right\rangle=\int_\mathbb{R} dx'\left(g'(x')f(x',t)+D(x',t)g''(x')\right)P(x';t|x_0;0)$

接下来对g(x)做分部积分。在以下三种通常使用的情况下边界项等于0：

周期性边界，两边界取值相等。
闭盒子，在盒子外粒子概率为0。
无穷空间，这时候概率密度可积保证了P(x,t)在无穷远处为0。

这时我们有

$\begin{align} \int_\mathbb{R} dx'\,\left(g'(x')f(x',t)+D(x',t)g''(x')\right)P(x';t|x_0;0)=&\int_\mathbb{R}dx'\, \left(-\frac{\partial f(x',t)P(x';t|x_0;0)}{\partial x'}+\frac{\partial^2D(x',t)P(x';t|x_0;0)}{\partial x'^2}\right)g(x') \\ =&\int_\mathbb{R} dx'\,g(x')\frac{\partial P(x';t|x_0,0)}{\partial t} \end{align}$

考虑到g(x)的任意性，这意味着

$\frac{\partial P(x;t|x_0,0)}{\partial t}=\frac{\partial^2D(x,t)P(x;t|x_0;0)}{\partial x^2}-\frac{\partial f(x,t)P(x;t|x_0;0)}{\partial x}$

这就是我们想要得到的对应于朗之万方程的主方程——（正向）Fokker-Planck方程。在高维情形，对于朗之万方程

$\frac{d \bm{X_t}}{dt}=\bm{f}(\bm{X}_t,t)+\bm{\sigma}(\bm{X}_t,t)\bm{\xi}(t)$

与之前完全一样的推导可以得到其FP方程是

$\frac{\partial P(\bm{x};t)}{\partial t}=\sum_{i,j,k=1}^d\nabla_i\nabla_j\frac{\sigma_{ik}\sigma_{jk}}{2}P-\sum_{i=1}^d \nabla_if_iP$

（特别（不严格）地可以取g(x)为微观的密度函数，注意x'是粒子的轨迹，x是空间坐标，

$\hat\rho(x,x',t)=\delta(x-x')$

因为 $P(x,t|x_0,0)=\langle\hat{\rho}(x-x'(t))\rangle_{x'}$ ，从这个记号可以更简单地推出FP方程。）

注意，这个FP方程是在伊藤积分的诠释下得到的。同样的一维朗之万方程按Stratonovich积分诠释，对应的FP方程应该是

$\frac{\partial P(x;t|x_0,0)}{\partial t}=\frac{\partial}{\partial x}\sqrt{D(x,t)}\left(\frac{\partial\sqrt{D(x,t)}P(x;t|x_0;0)}{\partial x}\right)-\frac{\partial f(x,t)P(x;t|x_0;0)}{\partial x}$

作为主方程，FP方程包含马尔可夫过程的全部信息，因此是与特定的随机过程——或者说物理——一一对应的。只有朗之万方程可以有不同的诠释。熟悉扩散方程的读者应该已经看出来了，FP方程实际上就是概率密度的扩散-漂移方程。定义扩散张量 $D_{ij}=\frac{1}{2}\sum_k\sigma_{ik}\sigma_{jk}$ ，我们可以把多维FP方程写成连续性方程的形式

$\frac{\partial P}{\partial t}=-\nabla\cdot\bm{J}$

$\bm{J}=\nabla(\bm{D}P)-\bm{f}P$

最后，与反向主方程类似地，我们也有反向Fokker-Planck方程

$-\frac{\partial P(x;t|x',t')}{\partial t'}=D(x',t')\frac{\partial^2P(x;t|x';t')}{\partial x'^2}+f(x',t')\frac{\partial P(x;t|x';t')}{\partial x'}$

其推导利用正向FP方程和马尔可夫过程的性质不难得出（此题老板上课时是留作作业的嗯）。

下面两节会稍微离题一下，做一点简明但较为深入的讨论。之后从第八节开始我们再回到引言中所说的平衡态统计的问题。

六、细致平衡条件

这里我们简单讨论一下FP方程描述的系统，其稳态的细致平衡条件。对于

$\dot{x}(t)=f(x)+\sqrt{2D(x)}\xi(t)$

细致平衡条件是说，联合概率密度 $P(x',t',x,t)=P(x,t',x',t)$ ，即在稳态的一定时间t'-t内，观察到从相空间中的一个点到另一个点的概率，等于其逆过程的概率。考虑到稳态条件，这也意味着

$P(x',t'|x,t)P_{SS}(x)=P(x,t'|x',t)P_{SS}(x')$

其中 $P_{SS}(x)$ 是稳态的概率分布。记FP方程的微分算子为

$\partial_t P(x;t)=-\hat{\mathcal{H}}(x)P(x;t)$

$\hat{\mathcal{H}}(x)=-\partial^2_xD+\partial_x f$

取t'=t+dt，非常接近于t，考虑到 $P(x';t|x;t)=\delta(x'-x)$ ，则

$(1-dt\hat{\mathcal{H}}(x))\delta(x-x')P_{SS}(x')=(1-dt\hat{\mathcal{H}}(x'))\delta(x'-x)P_{SS}(x)$

考虑到delta函数是偶函数，这意味着

$\hat{\mathcal{H}}(x)\delta(x-x')P_{SS}(x')=\hat{\mathcal{H}}(x')\delta(x-x')P_{SS}(x)$

在 $\mathbb{L}_2$ 内积 $\langle f,g\rangle=\int dx\, f(x)g(x)$ 下， $\hat{\mathcal{H}}(x)$ 的伴随算子满足 $\hat{\mathcal{H}}^\dagger(x')\delta(x-x')=\hat{\mathcal{H}}(x)\delta(x-x')$ ，容易验证

$\hat{\mathcal{H}}^\dagger(x)=-D\partial^2_x-f\partial_x$

（如果你还记得反向FP方程，这正是它对应的微分算子。）于是

$\hat{\mathcal{H}}(x)\delta(x-x')P_{SS}(x')=\hat{\mathcal{H}}(x)P_{SS}(x)\delta(x-x')$ （ $\langle P,\hat{\mathcal{H}}\delta\rangle=\langle \hat{\mathcal{H}}P,\delta\rangle$ ）

$\hat{\mathcal{H}}(x')\delta(x-x')P_{SS}(x)=P_{SS}(x)\hat{\mathcal{H}}(x')\delta(x-x')=P_{SS}(x)\hat{\mathcal{H}}^\dagger(x)\delta(x-x')$

细致平衡条件等价于

$\hat{\mathcal{H}}(x)P_{SS}(x)=P_{SS}(x)\hat{\mathcal{H}}^\dagger(x)$

如果f不含时，可积， $f=-\partial_x V$ ，而D是常数，那么容易验证 $P_{SS}\propto e^{-V/D}$ ，于是代入上式可以发现细致平衡条件是成立的。而同样的朗之万方程则具有时间反演对称性。一般情况下，细致平衡条件成立往往意味着朗之万方程具有时间反演对称性。

七、路径积分构造

这一节简单阐述一下如何将概率密度表示为路径积分的方式。首先考虑维那过程的概率密度函数。根据CK方程，

$P(x;t|x_0;t_0)=\int dx_1\,P(x;t|x_1;t_1)P(x_1;t_1|x_0;t_0)$

我们在其中取无数个分点

$P(x;t|x_0;t_0)=\int\cdots\int dx_1\cdots dx_n\,P(x;t|x_n;t_n)P(x_n;t_n|x_{n-1};t_{n-1})\cdots P(x_1;t_1|x_0;t_0)$

考虑到维那过程的增量服从正态分布，即 $P(x',t'|x,t)=e^{-(x'-x)^2/2(t'-t)}/\sqrt{2\pi(t'-t)}$ ，那么上式可以写成（记 $t_{n+1}=t,\ x_{n+1}=x$ ）

$P(x;t|x_0;t_0)=\int\cdots\int dx_1\cdots dx_n\,\prod_{k=1}^{n+1} \frac{1}{\sqrt{2\pi(t_k-t_{k-1})}}\exp\left(-\sum_{l=1}^{n+1}\frac{(x_l-x_{l-1})^2}{2(t_l-t_{l-1})}\right)$

在形式极限 $\max_k\{t_{k}-t_{k-1}\}\to0$ 下，前面的系数我们假装可以用一个归一化系数 $Z^{-1}$ 表示。而求和号内部的东西

$\sum_{l=1}^{n+1}\frac{(x_l-x_{l-1})^2}{2(t_l-t_{l-1})}=\frac{1}{2}\sum_{l=1}^{n+1}\left(\frac{x_l-x_{l-1}}{t_l-t_{l-1}}\right)^2(t_l-t_{l-1})=\frac{1}{2}\sum_l\dot{x}_l^2 \Delta t_l\sim\frac{1}{2}\int_0^t\dot{x}^2\, ds$