本文主要参考文献为 Robert Zwanzig 所著的 Nonequilibrium Statistical Mechanics.
1. 布朗运动与 Langevin 方程
1.1 Langevin 方程与涨落耗散定理
考虑一个球形粒子（半径为$~a~$，质量为$~m~$，坐标为$~x~$，速度为$~v~$）在液体介质中的一维运动，粒子运动的牛顿方程为
$$ m\frac{\mathrm d v}{\mathrm d t} = F_{\text{total}}(t) $$
其中$~F_{\text{total}}(t)~$是$~t~$时刻作用在粒子上的总的瞬时作用力。粒子与介质的作用是这个力的来源。在牛顿力学中，原则上这不是一个随机力。但想要分析其具体形式是很难做到的。实验告诉我们，在经典情形中，这个力主要是摩擦力，它的形式为$~-\zeta v~$，即与布朗粒子的速度成正比。由流体力学中的 Stokes 定律可知$~\zeta = 6\pi\eta a~$. 故式子变为
$$ m\frac{\mathrm d v}{\mathrm d t} \cong -\zeta v $$
其解为
$$ v(t) = e^{-\zeta t/m}v(0) $$
由解可以看出，长时间后粒子的速度应趋于零。但实际上处于热平衡的粒子的方均速度为
$$ \langle v^2\rangle_{\text{eq}} = \frac{kT}{m} $$
故我们对$~F_{\text{total}}(t)~$做出修正。我们给它加上称为随机力或涨落力的一项$~\delta F(t)~$，式子变为
$$ m\frac{\mathrm d v}{\mathrm d t} = -\zeta v + \delta F(t) $$
这就是布朗粒子的 Langevin 方程。假定随机力的一阶矩和二阶矩分别为
$$ \langle \delta F(t)\rangle = 0,\qquad \langle\delta F(t)\delta F(t')\rangle = 2B\delta(t-t') $$
常数$~B~$衡量涨落力的强度。$~\delta~$函数说明任意两个不同时刻的随机力是没有关联的。
Langevin 方程是一个一阶线性非齐次方程，它的解可以写为
$$ v(t) = e^{-\zeta t/m}v(0) + \int_0^t \mathrm d t'~\frac{\delta F(t')}{m}\exp\left[-\frac{\zeta(t-t')}{m}\right] $$
由一阶矩，二阶矩，以及方均速度的表达式可得
$$ B = \zeta kT $$
即
$$ \langle\delta F(t)\delta F(t')\rangle = 2\zeta kT\delta(t-t') $$
这个式子被称为涨落耗散定理。它说明涨落力的强度$~B~$与耗散的强度$~\zeta~$成正比，即涨落力越强，耗散越强。
1.2 时间关联函数
Langevin 方程可以用于计算各种时间关联函数。和平衡态统计力学相同，非平衡态统计力学也以系综理论为基础。但不同的时，处于非平衡态的系统并没有一个统一的能够计算各种量的配分函数。在非平衡态统计力学中，时间关联函数与配分函数、空间对关联函数有着相同的作用。
力学量$~A(t)~$的时间平均值为
$$ \lang A\rang = \frac{1}{\tau}\int_0^\tau\mathrm d t~A(t) $$
其涨落$~\delta A~$定义为
$$ \delta A(t) = A(t) - \langle A\rangle $$
不同时间的涨落可以是彼此相关的，在两个不同时刻的涨落乘积的时间平均值定义为涨落$~\delta A~$的时间关联函数（TCF）
$$ C(t) = \frac{1}{\tau}\int_0^\tau \mathrm d s~\delta A(s)\delta A(t+s) $$
我们常说的方均涨落，即是$~C(0)~$.
若所研究的系统满足各态历经假说，那么长时间的平均就等于系综平均，平衡态统计力学便有了用武之地。对于平衡态而言，我们可以通过计算两个涨落乘积的系综平均而不是其长时间平均来得到时间相关函数。在平衡系综中，不存在特殊的起始时间，$~C(t)~$的值仅取决于两次涨落之间的时间差。
对时间关联函数做傅里叶变换得到谱密度
$$ C_\omega = \int_{-\infty}^{\infty}\mathrm d t~e^{-\mathrm i \omega t}C(t) $$
谱密度通常可以由实验得到。对谱密度做反变换就得到了时间关联函数。
1.3 速度关联函数
时间关联函数的一个简单例子就是单粒子的速度关联函数。考虑一个一维的扩散，用$~C(x, t)~$表示粒子的浓度，该过程满足扩散方程
$$ \frac{\partial }{\partial t}C(x, t) = D\frac{\partial ^2}{\partial x^2}C(x, t) $$
计算粒子位置方均值的变化速率
$$ \begin{aligned} \frac{\partial }{\partial t}\langle x^2\rangle &= \int\mathrm d x~x^2\frac{\partial }{\partial t}C(x, t) \\ &= D\int\mathrm d x~x^2\frac{\partial ^2}{\partial x^2}C(x, t) \\ &=2D\int\mathrm d x~C(x, t) \\ &= 2D \end{aligned} $$
第二行到第三行由分部积分得到。最后假定浓度归一化到一。这就得到了著名的爱因斯坦一维扩散公式$~\langle x^2\rangle = 2Dt~$.
现在从统计的视角来看。粒子的位移为
$$ x(t) = \int_0^t\mathrm d s~v(s) $$
位移平方的系综平均为
$$ \begin{aligned} \langle x^2\rangle &= \left\langle\int_0^t\mathrm d s_1~v(s_1)\int_0^t\mathrm d s_2~v(s_2)\right\rangle \\ &= \int_0^t\mathrm d s_1\int_0^t\mathrm d s_2\langle v(s_1)v(s_2) \rangle \end{aligned} $$
求时间偏导得到
$$ \frac{\partial }{\partial t}\langle x^2\rangle = 2\int_0^t\mathrm d s~\langle v(t)v(s)\rangle $$
平衡态下的平均值应与时间轴的选取无关。该积分只与时间差$~t-s = u~$有关，故
$$ \begin{aligned} \frac{\partial }{\partial t}\langle x^2\rangle &= 2\int_0^t\mathrm d s~\langle v(t-s)v(0)\rangle \\ &= 2\int_0^t\mathrm d u~\langle v(u)v(0)\rangle \end{aligned} $$
对于极短的时间间隔，速度关联函数退化到零。故扩散方程只对长时间的情形有效。当$~t~$趋于无穷时，等式右边为$~2D~$，故有
$$ D = \int_0^\infty \mathrm d t~\langle v(t)v(0)\rangle $$
此式给出了扩散系数$~D~$与时间关联函数之间的关系。
对于三维情形，等式变为
$$ D = \frac{1}{3}\int_0^\infty\mathrm d t~\langle \boldsymbol V(t)\cdot\boldsymbol V(0) \rangle $$
其中$~\boldsymbol V~$为速度矢量。
1.4 关联函数与布朗运动
Langevin 方程和涨落耗散定理可以用来找到很多时间关联函数的表达式。
计算速度的平衡系综平均包括对噪声的平均和对初始速度的平均。噪声平均导致
$$ \langle v(t)\rangle_{\text{niose}} = e^{-\zeta t/m}v(0) $$
将$~v(t)~$与$~v(0)~$相乘再做平均得到
$$ \langle v(t)v(0)\rangle_{\text{eq}} = \frac{kT}{m}e^{-\zeta t/m} $$
此式只对$~t \ge 0~$有效，因为 Langevin 方程只对非负的时间成立。
我们希望速度关联函数是时间差的绝对值的函数，但要从 方程看这一点我们必须做长时间的平均。计算速度关联函数的长时间平均
$$ \langle v(t)v(t')\rangle_{\text{time}} = \frac{1}{\tau}\int_0^\tau\mathrm d s~v(t+s)v(t'+s) $$
$t~$时刻的瞬时速度由它的初值和对噪声的积分决定。我们假设初始时间是无穷远的过去，这样速度的初始值的贡献就衰减为零，瞬时速度就只由噪声决定了。得到
$$ v(t) = \int_0^\infty\mathrm d u~\frac{\delta F(t-u)}{m}e^{-\zeta u/m} $$
故速度关联函数是一个三重积分
$$ \begin{aligned} \langle v(t)v(t')\rangle_{\text{time}} &= \int_0^\infty\mathrm d u_1\int_0^\infty\mathrm d u_2~e^{-\zeta(u_1+u_2)/m}\frac{1}{\tau}\int_0^\tau\mathrm d s~\frac{1}{m^2}\delta F(t-u_1+s)\delta F(t-u_2+s) \\ &= \int_0^\infty\mathrm d u_1\int_0^\infty\mathrm d u_2~e^{-\zeta(u_1+u_2)/m}\frac{1}{\tau}\int_0^\tau\mathrm d s~\frac{2B}{m^2}\delta (t-u_1-t'+u_2) \\ &= \frac{kT}{m}e^{-\zeta|t-t'|/m} \end{aligned} $$
这里绝对值在计算过程中自动的出现了。
从 Langevin 方程的解入手，可以得到
$$ \langle x^2\rangle_{\text{eq}} = 2\frac{kT}{\zeta}\left[ t-\frac{m}{\zeta} + \frac{m}{\zeta}e^{-\zeta t/m} \right] $$
当$~t\to \infty~$时，有
$$ \langle x^2\rangle_{\text{eq}} \to 2\frac{kT}{\zeta}t $$
比较可得
$$ D = \frac{kT}{\zeta} $$
上式被称为 Stokes-Einstein 公式。
2. Fokker-Planck 方程
2.1 经典力学中的刘维尔方程
定义刘维尔算符
$$ L = \frac{\partial H}{\partial \boldsymbol p}\cdot\frac{\partial}{\partial \boldsymbol q} - \frac{\partial H}{\partial \boldsymbol q}\cdot\frac{\partial}{\partial \boldsymbol p} $$
则刘维尔方程可以写为
$$ \frac{\partial f}{\partial t} = -Lf $$
其中$~f~$是相空间的分布函数。
刘维尔方程有算符解
$$ f(\boldsymbol X, t) = e^{-tL}f(\boldsymbol X, 0) $$
下面给出$~L~$的一些性质。
考虑$~LAf~$在整个相空间上的积分，有
$$ \int_\omega\mathrm d \boldsymbol X~LAf = -\int_\omega\mathrm d \boldsymbol X~\frac{\mathrm d}{\mathrm d \boldsymbol X}\cdot(\boldsymbol V Af) = -\oint_{\partial \omega}\mathrm d \boldsymbol S~\cdot(\boldsymbol VAf) $$
其中$~\boldsymbol S~$是边界的法向量。系统有有限的能量，限制在位形空间中一块有限的区域中，故这样的积分变换是可以做的。
由于$~L~$含有微分运算，故$~L(Af) = (LA)f + A(Lf)~$. 另外，$~L~$在相空间中是反自伴的，即
$$ \int\mathrm d\boldsymbol X~A(Lf) = -\int\mathrm d \boldsymbol X~(LA)f $$
2.2 动力学变量
考虑由初始坐标$~\boldsymbol X~$描述的力学量$~A~$，它初始的变化速率为
$$ \begin{aligned} \left( \frac{\partial A}{\partial t} \right)_{t = 0} &= \frac{\partial A}{\partial \boldsymbol q}\cdot\left( \frac{\partial \boldsymbol q}{\partial t} \right)_{t = 0} + \frac{\partial A}{\partial \boldsymbol p}\cdot \left( \frac{\partial \boldsymbol p}{\partial t} \right)_{t = 0} \\ &= \left( \frac{\partial H}{\partial \boldsymbol p}\cdot\frac{\partial}{\partial \boldsymbol q} - \frac{\partial H}{\partial \boldsymbol q}\cdot\frac{\partial }{\partial p} \right)A \\ &= LA \end{aligned} $$
初始的变化率为$~LA~$. 可以证明，初始的二阶变化速率为$~L^2A~$. $~n~$阶初始变化速率为
$$ \left( \frac{\partial^n A}{\partial t^n} \right)_{t = 0} = L^n A $$
故由初始坐标$~\boldsymbol X~$描述的力学量$~A~$在$~t~$时刻时为
$$ \begin{aligned} A(\boldsymbol X, t) &= \sum_{j = 0}^\infty \frac{t^n}{n!}\left( \frac{\partial^n A}{\partial t^n} \right)_{t=0} \\ &= \sum_{j = 0}^\infty \frac{t^n}{n!}L^n A(\boldsymbol X) \\ &= e^{tL}A(\boldsymbol X) \end{aligned} $$
上式就是如下算符方程的解
$$ \frac{\partial }{\partial t}A(\boldsymbol X, t) = LA(\boldsymbol X, t), \qquad A(\boldsymbol X, 0) = A(\boldsymbol X) $$
这个方程就是动力学变量时间演化的刘维尔方程。
算符$~\exp(tL)~$给出了动力学变量在相空间中的演化轨迹，我们称此算符为传播子。
传播子有如下的一些性质。
它可以移到函数内
$$ e^{tL}A(\boldsymbol X) = A(e^{tL}\boldsymbol X) $$
它可以分配给两个相乘的函数
$$ e^{tL}A(\boldsymbol X)B(\boldsymbol X) = (e^{tL}A(\boldsymbol X))(e^{tL}B(\boldsymbol X)) $$
这是相空间轨迹互不相交的特性导致的。
变量$~A~$在$~t~$时刻的相空间平均为
$$ \langle A, t\rangle = \int\mathrm d \boldsymbol X~A(\boldsymbol X)f(\boldsymbol X, t) = \int\mathrm d \boldsymbol X~A(\boldsymbol X)e^{-tL}f(\boldsymbol X, 0) $$
这个平均还可以写为
$$ \langle A, t\rangle = \int\mathrm d \boldsymbol X~A(\boldsymbol X)f(\boldsymbol X, t) = \int\mathrm d \boldsymbol X~(e^{tL}A(\boldsymbol X))f(\boldsymbol X, 0) $$
这两个表达式是相等的，因为$~L~$在相空间是反自伴算子。
2.3 时间关联函数
时间关联函数可以写为
$$ \begin{aligned} C_{AB}(t) &= \int\mathrm d \boldsymbol X~A(\boldsymbol X, t)B(\boldsymbol X)f_{\text{eq}}(\boldsymbol X) \\ &=\int\mathrm d \boldsymbol X~(e^{tL}A(\boldsymbol X))B(\boldsymbol X)f_{\text{eq}}(\boldsymbol X) \\ &=\int\mathrm d \boldsymbol X~A(\boldsymbol X)(e^{-tL}B(\boldsymbol X)f_{\text{eq}}(\boldsymbol X)) \end{aligned} $$
指数算符可以用前面的性质分配给$~B~$和$~f~$，并且注意到平衡时的分布函数是不随时间变化而变化的，故有
$$ C_{AB}(t) = \int\mathrm d\boldsymbol X~A(t)B(\boldsymbol X, -t)f_{\text{eq}}(\boldsymbol X) = C_{AB}(-t) $$
如果$~A, B~$是同一个量，那么时间关联函数是时间反演不变的。
时间关联函数的导数是另一个时间关联函数
$$ \begin{aligned} \frac{\partial}{\partial t}C_{AB}(t) &= \int\mathrm d \boldsymbol X~\frac{\partial }{\partial t}A(\boldsymbol X, t)B(\boldsymbol X)f_{\text{eq}}(\boldsymbol X) \\ &= \int\mathrm d \boldsymbol X~(LA(\boldsymbol X, t))B(\boldsymbol X)f_{\text{eq}}(\boldsymbol X) \\ &= -\int\mathrm d \boldsymbol X~A(\boldsymbol X, t)(LB(\boldsymbol X))f_{\text{eq}}(\boldsymbol X) \end{aligned} $$
这是$~A~$与$~B~$的时间导数的时间关联函数。同理，有
$$ \frac{\partial^2}{\partial t^2}C_{AB}(t) = -\int\mathrm d \boldsymbol X~\dot{A}(\boldsymbol X, t)\dot{B}(\boldsymbol X)f_{\text{eq}}(\boldsymbol X) $$
例如，速度关联函数的二阶时间导数是负的随机力关联函数
$$ \frac{\partial^2}{\partial t^2}\langle v(t)v\rangle_{\text{eq}} = -\frac{1}{m^2}\langle F(t)F\rangle_{\text{eq}} $$
2.4 Fokker-Planck 方程
首先从 Langevin 方程入手。我们用符号$~\boldsymbol a~$表示一系列变量$~\{ a_1, a_2, \cdots \}~$. 我们要求过程是马氏性的，噪声是具有高斯分布的白噪声。则运动方程可以写为
$$ \frac{\mathrm d \boldsymbol a}{\mathrm d t} = \boldsymbol v(\boldsymbol a) + \boldsymbol F(t) $$
其中$~\boldsymbol v(\boldsymbol a)~$是变量$~\boldsymbol a~$的函数，噪声$~\boldsymbol F(t)~$是有高斯分布的，均值为零的白噪声，其二阶矩为
$$ \langle \boldsymbol F(t)\boldsymbol F(t')\rangle = 2\boldsymbol B\delta(t-t') $$
不同于 Langevin 方程，我们现在并不想找到这个方程的解，而是想看概率分布$~f(\boldsymbol a, t)~$的变化。更进一步的讲，我们想要知道概率分布的噪声平均。
首先认识到概率分布是一个守恒量
$$ \int\mathrm d \boldsymbol a~f(\boldsymbol a, t) = 1,\qquad\forall t $$
因此我们可以给出$~f~$的守恒律，即它的时间导数与流的散度是相制衡的。流可以写为速度乘上密度。我们令$~\boldsymbol a~$为空间坐标，则有守恒律
$$ \frac{\partial f}{\partial t} + \frac{\partial }{\partial \boldsymbol a}\cdot\left( \frac{\partial \boldsymbol a}{\partial t} f \right) = 0 $$
把$~\boldsymbol a~$的 Langevin 方程带入得到
$$ \frac{\partial f(\boldsymbol a, t)}{\partial t} = -\frac{\partial }{\partial \boldsymbol a}[\boldsymbol v(\boldsymbol a)f(\boldsymbol a, t) + \boldsymbol F(t)f(\boldsymbol a, t)] $$
这是一个随机微分方程。我们希望用它来求得$~f~$的噪声平均。就像处理刘维尔方程时一样，我们定义一个算符$~L~$。$~L~$作用于任何函数$~\Phi~$得到
$$ L\Phi \equiv \frac{\partial }{\partial \boldsymbol a}[\boldsymbol v(\boldsymbol a)\Phi] $$
故没有噪声时，方程可以写为
$$ \frac{\partial f}{\partial t} = -Lf $$
它有算符解
$$ f(\boldsymbol a, t) = e^{-tL}f(\boldsymbol a, 0) $$
现在加上噪声项
$$ \frac{\partial f}{\partial t} = -Lf - \frac{\partial }{\partial \boldsymbol a}[\boldsymbol F(t)f] $$
对时间积分得到
$$ f(\boldsymbol a, t) = e^{-tL}f(\boldsymbol a, 0) - \int_0^t\mathrm d s~e^{-(t-s)L}\frac{\partial}{\partial \boldsymbol a}[\boldsymbol F(s)f(\boldsymbol a, s)] $$
注意到$~f(\boldsymbol a, t)~$只依赖于$~t~$时刻之前的噪音$~\boldsymbol F(s)~$. 将上式回带到方程中得到
$$ \frac{\partial}{\partial t}f(\boldsymbol a, t) = -Lf(\boldsymbol a, t) - \frac{\partial}{\partial \boldsymbol a}[\boldsymbol F(t)f(\boldsymbol a, 0)] + \frac{\partial}{\partial \boldsymbol a}\boldsymbol F(t)\int_0^t\mathrm d s~e^{-(t-s)L}\frac{\partial}{\partial \boldsymbol a}[\boldsymbol F(s)f(\boldsymbol a, s)] $$
现在我们来做噪声平均。初始分布$~f(\boldsymbol a, 0)~$不受噪声平均的影响，所以单$~\boldsymbol F(t)~$项和初始分布函数的平均值为零。最后一项显式包含两个噪声因子$~\boldsymbol F(t)~$和$~\boldsymbol F(s)~$，以及$~f(\boldsymbol a, s)~$中隐式的噪声因子，但$~f(\boldsymbol a, s)~$中的噪声时间仅早于$~s~$。噪声为高斯函数，且相关函数是$~\delta~$函数。这意味着在做平均时，我们可以将第一个因子$~\boldsymbol F(t)~$与第二个因子$~\boldsymbol F(s)~$配对，或者与$~f(\boldsymbol a, s)~$中的一个隐式噪声因子配对。第一种情况下，我们得到$~\delta(t-s)~$；第二种情况下，我们得到$~\delta(t-s') ~$，其中$~s' < s~$. 但是第二种情况是不允许的，因为有$~t > s > s’~$的限制。因此，只需要对前两个噪声因子进行配对。平均值引入了因子$~\boldsymbol B~$，而$~\delta~$函数消去了算子$~e^{-(t-s)L}~$. 最终结果得到噪声平均分布函数$~\langle f(\boldsymbol a, t)\rangle~$的 Fokker-Planck 方程
$$ \frac{\partial}{\partial t}\langle f(\boldsymbol a, t)\rangle = -\frac{\partial}{\partial \boldsymbol a}[\boldsymbol v(\boldsymbol a)\langle f(\boldsymbol a, t)\rangle] + \frac{\partial }{\partial \boldsymbol a}\left[\boldsymbol B\frac{\partial}{\partial \boldsymbol a}\langle f(\boldsymbol a, t)\rangle\right] $$
右边第一项描述了没有噪声的情况，第二项则体现了噪声的影响。这里，$~\boldsymbol B~$可以是$~\boldsymbol a~$的任意函数。
这里给出的推导明确地依赖于两个假设。噪声是高斯噪声，它的相关函数是$~\delta~$函数。否则，对噪声的平均分解将不起作用。特别地，这种推导对于非马氏朗之万方程是无效的。
整个推到中我们并没有用到涨落耗散定理，也从没有说$~\langle f(\boldsymbol a, t)\rangle~$在长时间后一定会趋于平衡分布。
后面，我们会省略表示噪声平均的尖括号$~\langle~\rangle~$，我们只会讨论平均分布。
2.5 示例
考虑单粒子在势场$~U(x)~$中的布朗运动。Langevin 方程为
$$ \frac{\mathrm d x}{\mathrm d t} = \frac{p}{m},\qquad \frac{\mathrm d p}{\mathrm d t} = -U'(x) - \zeta \frac{p}{m} + F_p(t) $$
涨落耗散定理为
$$ \langle f_p(t)f_p(t') \rangle = 2\zeta kT\delta(t-t') $$
则 Fokker-Planck 方程方程中对应的变量为
$$ \begin{aligned} &\boldsymbol a = \begin{pmatrix} x \\ p \end{pmatrix},\qquad & \boldsymbol v(\boldsymbol a) = \begin{pmatrix} p/m \\ -U'(x) - \zeta p/m \end{pmatrix} \\ &\boldsymbol F(t) = \begin{pmatrix} 0 \\ F_p(t) \end{pmatrix},\qquad & \boldsymbol B = \begin{pmatrix} 0 & 0 \\ 0 & \zeta kT \end{pmatrix} \end{aligned} $$
相应的 Fokker-Planck 方程为
$$ \frac{\partial f}{\partial t} = -\frac{\partial}{\partial x}\frac{p}{m}f - \frac{\partial}{\partial p}\left[ -U'(x) - \frac{\zeta p}{m} \right]f + \zeta kT\frac{\partial^2 f}{\partial p^2} $$
注意到，如果没有噪音和摩擦，Fokker-Planck 方程就会退化为相空间中的刘维尔方程。哈密顿量为
$$ H = \frac{p^2}{2m} + U(x) $$
有噪声和摩擦存在时，方程的平衡解（令$~\partial f / \partial t = 0~$）为
$$ f_{\text{eq}}(x, p) = \frac{1}{Q}\exp\left[ -\frac{H(x, p)}{kT} \right] $$
其中$~Q~$为温度$~T~$时的配分函数
$$ Q = \iint\mathrm d x\mathrm d p~\exp\left[ -\frac{H(x, p)}{kT} \right] $$
假设系统的弛豫时间$~\tau = m/\zeta~$极其短，比任何与粒子在势场$~U(x)~$中运动有关的时间都要短很多。那么 Langevin 方程可以近似写为
$$ \frac{\mathrm d x}{\mathrm d t} = -\frac{1}{\zeta}U'(x) + \frac{1}{\zeta}F(t) $$
这个方程给出的 Fokker-Planck 方程我们通常称为斯莫鲁霍夫斯基方程
$$ \begin{aligned} \frac{\partial f}{\partial t} &= \frac{1}{\zeta}\frac{\partial}{\partial x}U'(x)f + \frac{kT}{\zeta}\frac{\partial^2 f}{\partial x^2} \\ &= D\frac{\partial}{\partial x}\exp\left[-\frac{U(x)}{kT}\right]\frac{\partial}{\partial x}\exp\left[\frac{U(x)}{kT}\right]f \end{aligned} $$
这个式子描述了外力驱动下的扩散，扩散系数为
$$ D = \frac{kT}{\zeta} $$
2.6 自伴算子
定义算符$~\mathcal D~$
$$ \mathcal D = -\frac{\partial}{\partial \boldsymbol a}\cdot \boldsymbol v(\boldsymbol a) + \frac{\partial }{\partial \boldsymbol a}\cdot \boldsymbol B\cdot \frac{\partial}{\partial \boldsymbol a} $$
则初始态为$~f(\boldsymbol a, t)~$的系统的时间演化满足 Fokker-Planck 方程
$$ \frac{\partial}{\partial t}f(\boldsymbol a, t) = \mathcal D f(\boldsymbol a, t) $$
算符$~\mathcal D~$的前一项就是我们之前定义过的刘维尔算符$~L~$，第二项则描述了噪声的平均作用。Fokker-Planck 方程有算符解
$$ f(\boldsymbol a, t) = e^{\mathcal D t}f(\boldsymbol a, 0) $$
这可以用于给出任意动力学变量$~\phi(\boldsymbol a)~$的平均值
$$ \langle \phi, t\rangle = \int\mathrm d \boldsymbol a~\phi(\boldsymbol a)f(\boldsymbol a, t) = \int\mathrm d \boldsymbol a~ \phi(\boldsymbol a)e^{\mathcal D t}f(\boldsymbol a, 0) $$
上式是我们把关注点放在概率分布的时间演化上，在$~t~$时刻求得的平均值。
求平均的第二种方法是用$~\mathcal D~$伴随算子，其定义为
$$ \begin{gathered} \int\mathrm d\boldsymbol a~\phi(\boldsymbol a)\mathcal D\Psi(\boldsymbol a) = \int\mathrm d \boldsymbol a~\Psi(\boldsymbol a)\mathcal D^\dagger\phi(\boldsymbol a) \\ \mathcal D^\dagger = \boldsymbol v(\boldsymbol a)\cdot\frac{\partial}{\partial \boldsymbol a} + \frac{\partial}{\partial \boldsymbol a}\cdot\boldsymbol B\cdot\frac{\partial}{\partial \boldsymbol a} \end{gathered} $$
则平均为
$$ \langle\phi,t \rangle = \int\mathrm d \boldsymbol a~f(\boldsymbol a, 0)e^{\mathcal D^\dagger t}\phi(\boldsymbol a) = \int\mathrm d \boldsymbol a~f(\boldsymbol a, 0)\phi(\boldsymbol a, t) $$
其中的含时变量$~\phi(\boldsymbol a, t)~$为
$$ \phi(\boldsymbol a, t) = e^{\mathcal D^\dagger t}\phi(\boldsymbol a) $$
这是我们把关注点放在动力学可观测量的时间演化上，对初始分布做的平均。$~\phi~$的运动方程为
$$ \frac{\partial}{\partial t}\phi(\boldsymbol a, t) = \left[ \boldsymbol v(\boldsymbol a)\cdot\frac{\partial}{\partial \boldsymbol a} + \frac{\partial}{\partial \boldsymbol a}\cdot\boldsymbol B\cdot\frac{\partial}{\partial\boldsymbol a} \right]\phi(\boldsymbol a, t),\qquad \phi(\boldsymbol a, 0) = \phi(\boldsymbol a) $$
量$~\phi(\boldsymbol a, t)~$的时间依赖性我们在实验上是看不到的，因为它是噪声平均前的情况。
2.7 格林函数
我们用指数形式的算符$~\exp(t\mathcal D)~$写出了 Fokker-Planck 方程的解。实际上，用格林函数可以$~G(\boldsymbol a, t|\boldsymbol a_0)~$可以写出更明确的解。
$$ f(\boldsymbol a, t) = e^{t\mathcal D}f(\boldsymbol a, 0) = \int\mathrm d\boldsymbol a~G(\boldsymbol a, t|\boldsymbol a_0)f(\boldsymbol a_0, 0) $$
这个格林函数是有着特殊初始条件的 Fokker-Planck 方程的解
$$ G(\boldsymbol a, t=0|a_0) = \delta(\boldsymbol a - \boldsymbol a_0) $$
如果流函数$~\boldsymbol v(\boldsymbol a)~$是$~\boldsymbol a~$的线性函数，即$~\boldsymbol v(\boldsymbol a) = \Theta\cdot\boldsymbol a~$，那么格林函数可以用下面的操作很容易的找到。首先对$~G~$做傅里叶变换
$$ \gamma(\xi, t|\boldsymbol a_0) = \int\mathrm d \boldsymbol a~e^{\mathrm i\xi\cdot \boldsymbol a}G(\boldsymbol a, t|\boldsymbol a_0) $$
然后用分部积分，Fokker-Planck 方程变成一阶偏微分方程
$$ \frac{\partial}{\partial t}\gamma = \xi\cdot\Theta\cdot\frac{\partial}{\partial \xi}\gamma - \xi\cdot\boldsymbol B\cdot\xi\gamma $$
则$~\gamma~$的对数满足方程
$$ \frac{\partial}{\partial t}\ln\gamma = \xi\cdot\Theta\cdot\frac{\partial}{\partial \xi}\ln\gamma - \xi\cdot\boldsymbol B\cdot\xi $$
让$~\ln\gamma~$对$~\xi~$做泰勒展开，保留到二阶项，即
$$ \ln\gamma = \mathrm i\alpha(t)\cdot\xi - \frac{1}{2}\xi\cdot\sigma(t)\cdot\xi $$
两个时间依赖的系数$~\alpha~$和$~\sigma~$分别满足方程
$$ \frac{\partial}{\partial t}\alpha = \Theta\cdot\alpha,\qquad\frac{\partial}{\partial t}\sigma = \Theta\cdot\sigma + \sigma\cdot\Theta^\dagger + 2\boldsymbol B $$
由于$~\gamma~$的初值为$~\exp(\mathrm i\xi\cdot\boldsymbol a_0)~$，故系数的初值为
$$ \alpha(0) = \boldsymbol a_0,\qquad \sigma(0) = 0 $$
可以解得
$$ \alpha(t) = \langle\boldsymbol a; t\rangle,\qquad\sigma(t) = \langle(\boldsymbol a - \langle\boldsymbol a; t\rangle)(\boldsymbol a - \langle\boldsymbol a; t\rangle);t\rangle $$
再做傅里叶逆变换得到
$$ G(\boldsymbol a, t|\boldsymbol a_0) = \frac{1}{\sqrt{\det(2\pi\sigma(t))}}\exp\left[ -\frac{1}{2}\left(\boldsymbol a - e^{t\Theta}\cdot\boldsymbol a_0\right)\cdot\sigma(t)^{-1}\cdot\left(\boldsymbol a - e^{t\Theta}\cdot\boldsymbol a_0\right) \right] $$
前面的系数是归一化因子。
3. 由过阻尼 Langevin 方程推导 Fokker-Planck 方程
本节给出一种较为简单且直观的方法从过阻尼 Langevin 方程推导 Fokker-Planck 方程。
考虑过阻尼 Langevin 方程（即忽略速度的变化）
$$ \dot{\boldsymbol r} = -\eta^{-1}\boldsymbol \nabla V + \sqrt{\frac{2k_B T}{\eta}}\boldsymbol \xi(t). $$
考虑$~t~$时刻在位置$~\boldsymbol x~$找到粒子的概率
$$ \rho(\boldsymbol x, t) = \int\mathrm d \boldsymbol x \rho(\boldsymbol x, t)\delta^3(\boldsymbol x - \boldsymbol r(t)) = \langle \delta^3(\boldsymbol x - \boldsymbol r(t)) \rangle, $$
其中$~\boldsymbol x~$是一个确定的位置，$~\boldsymbol r~$具有随机性。
记 Langevin 方程为
$$ \frac{\mathrm d \boldsymbol r(t)}{\mathrm d t} = \boldsymbol v(\boldsymbol r(t)) + \boldsymbol u(t), $$
其中
$$ \boldsymbol v = -\eta^{-1}\boldsymbol \nabla V, \quad \boldsymbol u = \sqrt{\frac{2k_B T}{\eta}}\boldsymbol \xi(t). $$
考虑一段极短的时间$~\Delta t~$，位置变化量为
$$ \Delta \boldsymbol r(t) = \boldsymbol r(t + \Delta t) - \boldsymbol r(t) = \int_{t}^{t + \Delta t}\mathrm d \tau \boldsymbol v(\boldsymbol r(\tau)) + \int_t^{t+\Delta t}\mathrm d \tau \boldsymbol u(\tau). $$
对$~t + \Delta t~$时刻的$~\delta~$函数做泰勒展开，保留至二阶项
$$ \begin{aligned} \delta^3(\boldsymbol x - \boldsymbol r(t + \Delta t)) ={}& \delta^3(\boldsymbol x - \boldsymbol r(t)) - \Delta \boldsymbol r_i(t)\partial_i\delta^3(\boldsymbol x - \boldsymbol r(t)) \\ {}&+ \frac{1}{2}\Delta\boldsymbol r_i(t)\Delta\boldsymbol r_j(t)\partial_i\partial_j\delta^3(\boldsymbol x - \boldsymbol r(t)). \end{aligned} $$
其中$~\partial_i \equiv \partial / \partial x_i~$。
概率密度的变化量为
$$ \begin{aligned} \rho(\boldsymbol x, t + \Delta t) - \rho(\boldsymbol x, t) &= \langle \delta^3(\boldsymbol x, \boldsymbol r(t + \Delta t)) \rangle - \langle \delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle \\ &= -\partial_i\langle \Delta r_i(t)\delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle + \frac{1}{2}\partial_i\partial_j\langle \Delta r_i(t)\Delta r_j(t)\delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle \end{aligned} $$
上式左边做泰勒展开得到
$$ \rho(\boldsymbol x, t + \Delta t) - \rho(\boldsymbol x, t) \approx \Delta t\partial_t \rho(\boldsymbol x, t). $$
现在计算等式右边的两项。第一项为
$$ \begin{aligned} \langle \Delta r_i(t)\delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle &= \langle\Delta r_i(t)\rangle\langle \delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle \\ &= \langle \Delta r_i(t)\rangle \rho(\boldsymbol x, t). \end{aligned} $$
其中
$$ \langle \Delta r_i(t)\rangle = \int_t^{t+\Delta t}\mathrm d \tau\left( -\frac{1}{\eta}\langle \boldsymbol \nabla V \rangle \right) + \int_t^{t + \Delta t}\mathrm d \tau \sqrt{\frac{2k_B T}{\eta}}\langle\xi(t)\rangle = -\frac{\Delta t}{\eta}\boldsymbol \nabla V. $$
第二项为
$$ \langle \Delta r_i(t)\Delta r_j(t)\delta^3(\boldsymbol x, \boldsymbol r(t)) \rangle = \langle\Delta r_i(t)\Delta r_j(t)\rangle \langle \delta^3(\boldsymbol x, \boldsymbol r(t))\rangle. $$
其中
$$ \begin{aligned} \langle\Delta r_i(t)\Delta r_j(t)\rangle &= \int_t^{t + \Delta t}\mathrm d t_1\int_t^{t+\Delta t}\mathrm d t_2(\boldsymbol v_1\boldsymbol v_2 + \boldsymbol v_1\boldsymbol u_2 + \boldsymbol v_2\boldsymbol u_1 + \boldsymbol u_1\boldsymbol u_2) \\ &= \frac{\Delta t^2}{\eta^2}(\boldsymbol \nabla V)^2 + \frac{2k_B T}{\eta}\Delta t\delta_{ij}. \end{aligned} $$
最后令$~\Delta t \to 0~$，得到过阻尼 Langevin 方程对应的 Fokker-Planck 方程
$$ \partial_t \rho(\boldsymbol x, t) = \nabla\left[ \frac{1}{\eta}\boldsymbol \nabla V \rho(\boldsymbol x, t) \right] + \frac{k_B T}{\eta}\nabla^2 \rho(\boldsymbol x, t). $$
令概率流
$$ J(\boldsymbol x) = -\frac{\rho}{\eta}\boldsymbol \nabla V + \frac{k_B T}{\eta}\boldsymbol \nabla \rho, $$
则方程可以简写为
$$ \frac{\partial\rho}{\partial t} = -\nabla J. $$
