软物质的统计力学（四）

八、回到平衡态统计

有了概率密度，有了主方程，现在我们终于可以开始回答引言中提到的问题了——液体中悬浮的粒子这种通常用朗之万方程描述的系统还能不能使用正则系综？为此我们回想第二节中我们写下的在外场中、无相互作用的粒子的朗之万方程

$\begin{align} \frac{d\bm{q}}{dt}=&\frac{\bm{p}}{M}\\ \frac{d\bm{p}}{dt}=&-\nabla V(\bm{q})-\frac{\mu}{M}\bm{p}+\sqrt{2\mu k_BT}\bm{\xi}(t) \end{align}$

根据第五节给出的结果，我们现在可以轻易地写出它的概率密度 $P(\bm{q},\bm{p};t)$ 满足的FP方程

$\partial_tP(\bm{q},\bm{p};t)=-\frac{1}{M}\nabla_\bm{q}\cdot(\bm{p}P)+\nabla_\bm{p}\cdot\left(\frac{\mu}{M}\bm{p}P+\nabla_\bm{q}V(\bm{q})P\right)+\nabla^2_\bm{p}(\mu k_BTP)$

要验证它是否满足正则系综其实很简单。考虑这个系统的稳态， $\partial_tP=0$ ，我们将正则系综 $P=Z^{-1}e^{-H(\bm{q},\bm{p})/k_BT}$ 代入验证一下是不是方程的解就可以了。其中哈密顿量H就是正常的哈密顿量

$H(\bm{q},\bm{p})=\frac{p^2}{2M}+V(\bm{q})$

考虑到Z只是一个常数，我们不用管它，于是

$-\frac{1}{M}\nabla_\bm{q}\cdot(\bm{p}P)=-\frac{\bm{p}}{M}\cdot\nabla_\bm{q}P=\frac{\bm{p}\cdot\nabla_\bm{q}H}{ZMk_BT}e^{-H/k_BT}=\frac{\bm{p}\cdot\nabla_\bm{q}V}{ZMk_BT}e^{-H/k_BT}$

$\nabla_\bm{p}\cdot\frac{\mu}{M}\bm{p}P=\frac{d\mu}{M}P+\frac{\mu\bm{p}}{M}\cdot\nabla_\bm{p}P=\frac{d\mu}{ZM}e^{-H/k_BT}-\frac{\mu p^2}{ZM^2k_BT}e^{-H/k_BT}$

$\nabla_\bm{p}\cdot(\nabla_\bm{q}V(\bm{q})P)=-\frac{\bm{p}\cdot\nabla_\bm{q}V}{ZMk_BT}e^{-H/k_BT}$

$\nabla^2_\bm{p}(\mu k_BTP)=-\mu k_BT\nabla_\bm{p}\cdot(\frac{\bm{p}}{ZMk_B T}e^{-H/k_BT})=-\frac{d\mu}{ZM}e^{-H/k_BT}+\frac{\mu p^2}{ZM^2k_BT}e^{-H/k_BT}$

把它们全部加起来，你可以看到每一项都被精确地抵消。也就是说， $P=Z^{-1}e^{-H(\bm{q},\bm{p})/k_BT}$ 的确就是系统的稳态解，而且这个稳态解不与阻尼系数 $\mu$ 有关——也就是说在一个过阻尼系统中这个结论依然成立，只是这种情况下动量的部分就直接被积分掉了。这符合我们对阻尼系数的概念——它是一个过程参量，不应该参与影响稳态（当然，更普遍的情况下这个结论不成立）。

那么，既然概率分布函数依然保持了正则系综的形式，那么基于概率分布函数定义的一切物理量——熵、自由能等等，以及基于这些函数的所有分析，都回到了仿佛没有流体的平衡态统计之中，尽管我们的运动方程因为阻尼项和噪声项的存在，并不是一个经典的哈密顿系统！之所以可以回到平衡态，大概就是因为涨落耗散定理使得系统的涨落和耗散有紧密的联系，使得它们对概率密度的贡献恰好抵消。于是我们看到，在没有相互作用的情况下，朗之万方程描述的粒子的统计依然是满足正则系综的。下一节我们则要考虑多粒子相互作用的情况。

九、平均场近似

现在我们考虑相互作用的N个全同粒子，第i个粒子的朗之万方程为

$\begin{align} \frac{d\bm{q}_i}{dt}=&\frac{\bm{p}_i}{m}\\ \frac{d\bm{p}_i}{dt}=&-\nabla V_{\mathrm{ext}}(\bm{q}_i)-\sum_{j\neq i}\nabla_{\bm{q}_i} V(\bm{q}_i-\bm{q}_j)-\frac{\mu}{m}\bm{p}_i+\sqrt{2\mu k_BT}\bm{\xi}_i(t) \end{align}$

其中 $V_{\mathrm{ext}}$ 是外场的势， $V$ 是相互作用势。我们当然可以写出N粒子概率密度函数 $P(\{\bm{q}_i,\bm{p}_i\};t)$ 满足的FP方程，但考虑到d维空间中P有2dN+1个变量，不用写出来也知道这个方程会非常长，而且因为相互作用造成的耦合，它不能表示为单粒子概率密度的积，我们很难处理这个大方程。我们转而考虑粒子的微观密度

$\hat{\rho}(\bm{q},\bm{p};t)=\sum_{i=1}^N\delta(\bm{q}-\bm{q_i}(t))\delta(\bm{p}-\bm{p_i}(t))$

而实验中可以观测到的，是这个函数的系综平均，我们称为（宏观）密度（注意，第二行积分中的变量全是相空间中的坐标，不再表示粒子的轨迹）

$\begin{align} \rho(\bm{q},\bm{p};t)=&\left\langle\sum_{i=1}^N\delta(\bm{q}-\bm{q_i}(t))\delta(\bm{p}-\bm{p_i}(t))\right\rangle \\ =&\int\cdots\int d\bm{q}_1\cdots d\bm{q}_Nd\bm{p}_1\cdots d\bm{p}_N P(\{\bm{q}_i,\bm{p}_i\};t)\sum_{i=1}^N\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i}) \\ =&\sum_{i=1}^N\int\cdots\int d\bm{q}_1\cdots d\bm{q}_{i-1}d\bm{q}_{i+1}\cdots d\bm{q}_Nd\bm{p}_1\cdots d\bm{p}_{i-1}d\bm{p}_{i+1}\cdots d\bm{p}_N P(\{\bm{q}_j,\bm{p}_j|\bm{q}_i=\bm{q},\bm{p}_i=\bm{p}\};t) \end{align}$

这个函数展开后看起来非常吓人，不过我们不用算它。我们可以利用伊藤积分计算密度随时间的演化。首先有（等式第二项利用了 $\partial\delta(\bm{q}-\bm{q}_j)/\partial\bm{q}_i=0$ 如果 $i\neq j$ ） $\begin{align} \frac{d\hat{\rho}}{dt}=&\sum_i\left(\frac{\partial\hat\rho}{\partial\bm{q}_i}\cdot\dot{\bm{q}}_i+\frac{\partial\hat\rho}{\partial\bm{p}_i}\cdot\dot{\bm{p}}_i\right)+\mu k_BT\sum_{i,j}\frac{\partial^2\hat{\rho}}{\partial\bm{p}_i\partial\bm{p}_j} \\ =&\sum_i\left(\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{q}_i}\cdot\dot{\bm{q}}_i+\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{p}_i}\cdot\dot{\bm{p}}_i+\mu k_BT\frac{\partial^2\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{p}_i^2}\right) \\ =&-\sum_i\left(\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{q}}\cdot\dot{\bm{q}}_i+\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{p}}\cdot\dot{\bm{p}}_i+\mu k_BT\frac{\partial^2\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{p}^2}\right) \end{align}$

代入朗之万方程，我们得到 $\begin{align} \frac{d\hat{\rho}}{dt}=&-\sum_i\left(\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{q}}\cdot\frac{\bm{p}_i}{m}\right. \\ &\left.+\frac{\partial\delta(\bm{q}-\bm{q_i})\delta(\bm{p}-\bm{p_i})}{\partial\bm{p}}\cdot\left(-\nabla V_{\mathrm{ext}}(\bm{q}_i)-\sum_{j\neq i}\nabla_{\bm{q}_i} V(\bm{q}_i-\bm{q}_j)-\frac{\mu}{m}\bm{p}_i+\sqrt{2\mu k_BT}\bm{\xi}_i(t)\right)\right)+\mu k_BT\frac{\partial^2\hat{\rho}}{\partial\bm{p}^2} \end{align}$

然后我们对方程做系综平均（实际上就是前面那个2dN重积分）。注意到 $\int \delta(x-x')f(x')\,dx'=\int \delta(x-x')f(x)\,dx'$ ， $\int f(x')\partial_x\delta(x-x')\,dx'=\int \partial_x(\delta(x-x')f(x))\,dx'$ ，系综平均中的积分中不含 $\bm{p}$ 和 $\bm{q}$ ，以及伊藤积分保证了噪声项的平均是0，再利用

$\sum_{j|j\neq i}\nabla_{\bm{q}_i}V(\bm{q}_i-\bm{q}_j)=\int d\bm{q}'\, \nabla_{\bm{q}_i}V(\bm{q}_i-\bm{q}')\sum_{j|j\neq i}\delta(\bm{q}'-\bm{q}_j)$

我们得到

$\begin{align} \frac{d\rho}{dt}=&-\frac{\partial}{\partial\bm{q}}\cdot\frac{\rho\bm{p}}{m}-\frac{\partial}{\partial\bm{p}}\cdot\left(-\nabla V_{\mathrm{ext}}(\bm{q})\rho-\frac{\mu}{m}\bm{p}\rho\right. \\ &\left.-\int d\bm{q}'d\bm{p}'\,\nabla_{\bm{q}} V(\bm{q}-\bm{q}')\left\langle\sum_{i,j|i\neq j}\delta(\bm{q}-\bm{q}_i)\delta(\bm{p}-\bm{p}_i)\delta(\bm{q}'-\bm{q}_j)\delta(\bm{p}'-\bm{p}_j))\right\rangle\right)+\mu k_BT\frac{\partial^2\rho}{\partial\bm{p}^2} \\ =&-\frac{\partial}{\partial\bm{q}}\cdot\frac{\rho\bm{p}}{m}+\frac{\partial}{\partial\bm{p}}\cdot\left(\nabla V_{\mathrm{ext}}(\bm{q})\rho+\int d\bm{q}'d\bm{p}'\,\nabla_{\bm{q}} V(\bm{q}-\bm{q}')\rho^{(2)}(\bm{q},\bm{p},\bm{q'},\bm{p'};t)+\frac{\mu}{m}\bm{p}\rho\right)+\mu k_BT\frac{\partial^2\rho}{\partial\bm{p}^2} \end{align}$

这里我们引入了两点粒子密度

$\rho^{(2)}(\bm{q},\bm{p},\bm{q'},\bm{p'};t)=\left\langle\sum_{i,j|i\neq j}\delta(\bm{q}-\bm{q}_i)\delta(\bm{p}-\bm{p}_i)\delta(\bm{q}'-\bm{q}_j)\delta(\bm{p}'-\bm{p}_j))\right\rangle$

为在相空间内两个点看到的粒子密度。虽然这个方程看起来比原来简单了一点，但因为我们不知道两点粒子密度的演化是什么，所以它并不是一个封闭的方程。诚然我们可以继续利用主方程将两点粒子密度的时间演化写下来，它肯定会与三点粒子密度相关。然后我们一步步走下去，会得到无穷多个微分方程——然后没有办法解了。所以，我们又到了要做近似的时候了。这个近似就是平均场近似

$\rho^{(2)}(\bm{q},\bm{p},\bm{q}',\bm{p}';t)=\rho(\bm{q},\bm{p};t)\rho(\bm{q}',\bm{p}';t)$

即假设相空间上两点的概率密度不相关。于是我们终于得到了一个闭合的方程来描述粒子的密度

$\frac{d\rho}{dt}=-\frac{\partial}{\partial\bm{q}}\cdot\frac{\rho\bm{p}}{m}+\frac{\partial}{\partial\bm{p}}\cdot\left(\nabla V_{\mathrm{ext}}(\bm{q})+\int d\bm{q}'d\bm{p}'\,\nabla_{\bm{q}} V(\bm{q}-\bm{q}')\rho(\bm{q'},\bm{p'};t)+\frac{\mu}{m}\bm{p}\right)\rho+\mu k_BT\frac{\partial^2\rho}{\partial\bm{p}^2}$

当然，这个方程既不是线性的，也不是局域的。但是，我们又可以猜方程的稳态解或许满足正则系综。现在我们取哈密顿量

$H(\bm{q},\bm{p})=\frac{p^2}{2m}+V_{\mathrm{ext}}(\bm{q})+\int d\bm{q}'d\bm{p}'\,V(\bm{q}-\bm{q'})\rho(\bm{q}',\bm{p'})$

然后我们猜 $\rho(\bm{q},\bm{p})=Z^{-1}e^{-H(\bm{q},\bm{p})/k_BT}$ 是系统的稳态解（当然，现在正则系综只给出了 $\rho$ 满足的一个积分方程）。接下来的计算和上一节完全一样，只是更长一点，我们很容易验证正则分布给出的的确就是方程的稳态解。于是，我们又回到了和上节一样的结论，在平均场近似下，有相互作用朗之万方程粒子系统也可以用平衡态正则系综来处理，且液体阻尼不影响平衡态结果。

前面的近似做得太快，以至于我们现在需要仔细思考一下刚才发生了什么事情。这个数学上没有道理的近似，在物理上对应着什么样的考虑？首先相空间内的粒子密度要足够高，而且势函数足够平缓，以至于我们去掉 $i\neq j$ 的限制——即加上了一个粒子对自己的“相互作用”——而引入的误差小到可以被忽略。其次，相空间上两点的密度不相关，意味着这两点之间几乎没有相互作用或关联。换句话说，空间两点之间的距离远大于系统的特征关联尺度 $\eta$ （粗略地说， $\langle\rho(x)\rho(x+l)\rangle\sim e^{-l/\eta}$ ）。这个条件在两种情况下被破坏：长程相互作用和相变临界点。如果我们选取的描述系统的尺度不合适，使得相互作用距离相对于空间尺度太大，那么显然平均场假设不能随意使用。而在相变临界点，众所周知的是，尽管相互作用是短程的，系统内的关联长度却是发散的，这时候空间任意两点的概率辐都是相关的，平均场近似也就失效了。事实也确实表明，平均场理论对于临界现象是不起作用的。而在两种情况下还能否回到平衡态正则系综，也是尚不清楚的。

在平均场近似之后，我们复杂的相互作用势变成了一个类似外场一样的东西

$V_{\mathrm{MF}}(\bm{q})=\int d\bm{q}'d\bm{p}'\,V(\bm{q}-\bm{q'})\rho(\bm{q}',\bm{p'})$

粒子感受到的平均场，是周围粒子在平均的意义下给自己的贡献之和。在抹平了周围粒子的贡献之后，多粒子系统的密度演化方程重新变得好像是单粒子一样了。而这种用更大的尺度来平滑点粒子场，进而用平滑的函数来描述粒子系统的方式，我们称为粗粒化。

接下来我们就可以像平时一样，用热力学来研究粒子系统，并计算相变发生时的相平衡条件。然而遗憾的是，通常用作描述分子间相互作用的Lennard-Jones势，性质并没有好到可以让我们直接使用前述的平均场近似，因为加入粒子与自己的相互作用必然得到发散的结果。这意味着我们可能需要更复杂的技巧来处理这类问题，而之后相关的问题也将缺少解析解。而究竟如何处理以LJ势相互作用的粒子，这是我现在还不清楚的。

来源：知乎 www.zhihu.com

作者：赵永峰

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载

闪新闻

软物质的统计力学（四）

目录

八、回到平衡态统计

九、平均场近似