经典网文分享

大坑的最后一篇，终于写完了。其实这几篇只包括了最基础的内容，感觉自己终于接近了大神们本科时的水平。虽然其实并没讲到非平衡态统计，也有很多理解不到位的地方，不过但愿这个系列能对想要做非平衡态统计的读者（比如我）有所帮助。

================================

十、密度泛函理论

现在我们只考虑平均场近似下有相互作用的粒子系统，在拓展一些概念的同时验证一下热力学定律。这时之前定义的粒子密度 $\rho(\bm{r},\bm{p})$ 其实某种意义下可以认为正比于在相空间 $(\bm{r},\bm{p})$ 处找到粒子的概率。那么我们可以仿照一般统计力学中的定义，定义熵

$s=-k_b\rho\ln\rho$

那么（亥姆霍兹）自由能密度 $f$ 依然可以写成

$f=h-Ts$

其中h是能量密度。

在实践中我们常常只关心粒子在位型空间中的密度，因此我们可以将动量的部分积分掉。这其实相当于考虑一个过阻尼的系统——如果粒子感受到的阻尼非常强，我们总可以选取一个足够长的时间尺度，使得粒子的动量变化可以被忽略 $\dot{\bm{p}}_i=0$ 。而我们之前也指出过系统的稳态性质实际上与阻尼系数无关，所以过阻尼系统的平衡态热力学应该是与一般系统一样的。在这种情况下，我们可以忽略粒子运动的动能部分。其朗之万方程简化为

$\frac{d\bm{r}_i}{dt}=-\frac{1}{\mu}\nabla V_{\mathrm{ext}}(\bm{r}_i)-\frac{1}{\mu}\sum_{j\neq i}\nabla_{\bm{r}_i} V(\bm{r}_i-\bm{r}_j)+\sqrt{\frac{2 k_BT}{\mu}}\bm{\xi}_i(t)$

而平均场近似下，粒子数密度 $\rho(\bm{r})$ 满足的运动方程简化为

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

既然忽略了动量部分，那么系统的能量密度可以写为所有势能的贡献之和（注意任意两个点之间的势能在积分时被计算了两次，所以要除以2）

$h(\bm{r})=V_{\mathrm{ext}}(\bm{r})\rho(\bm{r})+\frac{1}{2}\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'$

那么我们有自由能密度

$f(\bm{r})=V_{\mathrm{ext}}(\bm{r})\rho(\bm{r})+\frac{1}{2}\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'+k_BT\rho(\bm{r})\ln\rho(\bm{r})$

而它对整个空间的积分，就是整个系统的自由能（文献中更常见的做法是把外场拿出来放在自由能外面，这样我们可以把粒子间相互作用和外场分开来处理。这里为了行文方便没有这样做）

$\mathcal{F}[\rho]=\int_{\mathbb{R}^d}\left(V_{\mathrm{ext}}(\bm{r})\rho(\bm{r})+\frac{1}{2}\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'+k_BT\rho(\bm{r})\ln\rho(\bm{r})\right)\,d^d\bm{r}$

对于一个给定的密度函数 $\rho(\bm{r})$ ，我们计算出的整个系统的自由能是一个实数，也就是说，系统的自由能称为了密度函数的泛函。那么平衡态时密度还是否还对应于自由能泛函的极小？首先我们可以把 $\rho(\bm{r})$ 满足的运动方程写成

$\frac{\partial\rho}{\partial t}=\nabla\cdot\left(\frac{\rho}{\mu}\nabla\frac{\delta\mathcal{F}[\rho]}{\delta\rho}\right)$

也就是说，我们发现粒子的流

$\bm{J}=-\frac{\rho}{\mu}\nabla\frac{\delta\mathcal{F}[\rho]}{\delta\rho}$

平衡态时我们要求系统没有宏观流动 $\bm{J}=0$ （注：没有宏观流动主要指位型空间没有宏观流动，而不是相空间没有宏观流动。而且，实际上流的定义有一个"规范"自由度——所有可以表示为 $\bm{J}'=\bm{J}+\nabla\times\bm{A}$ 的流实际上是等价的。恰当地选取A也可以让平衡态时的流在相空间没有流动。）那么这意味着

$\nu=\frac{\delta\mathcal{F}[\rho]}{\delta\rho}$

是一个与空间坐标无关的常数。如果回想起化学势是等温等体条件下单位粒子的自由能，这里我们可以称 $\nu(\bm{r})$ 为自由能密度。

那么自由能密度等于常数是否就意味着这时候的自由能取道极值呢？自由能——或者说势能，其实有零点选取的自由。这意味着我们可以随意给 $\mathcal{F}[\rho]$ 加上一个常数。而实际上因为粒子数守恒，我们对粒子的密度有一个额外的限制条件 $\int_{\mathbb{R}^d}\rho\, d^d\bm{r}=N$ 是常数。也就是说，自由能

$\mathcal{F}'[\rho]=\mathcal{F}[\rho]+\int_{\mathbb{R}^d} A\rho(\bm{r})\, d^d\bm{r}$

与原自由能是等价的。但是，

$\frac{\delta\mathcal{F}'[\rho]}{\delta\rho}=\frac{\delta\mathcal{F}[\rho]}{\delta\rho}+A$

因此我们总可以恰当地选取A使得化学势等于0，因此化学势为常数就意味着自由能取到了极值。换句话说，化学势的绝对值也没有意义。

如果我们考察的空间尺度足够大，势能和密度的变化足够平滑，我们还可以将相互作用势展开，即

$\begin{align} \int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'=&\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}')\rho(\bm{r}-\bm{r}')\,d^d\bm{r}' \\ =&\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}')\left(\rho(\bm{r})-\bm{r'}\cdot\nabla\rho(\bm{r})+\frac{1}{2}\sum_{\alpha,\beta}r'_\alpha r'_\beta\partial_\alpha\partial_\beta\rho+o(\nabla^3\rho)\right)\,d^d\bm{r}' \end{align}$

如果 $V(\bm{r})=V(r)$ 是一个球对称的势函数，那么对称性使得 $\int_\mathbb{R} r_\alpha V(r)\,dr_\alpha=0$ ，所有的奇数次项都是0。于是，我们记

$\gamma=\frac{1}{2}\int_{\mathbb{R}^d}V(r)\,d^d\bm{r}$

$\kappa=-\frac{1}{2}\int_{\mathbb{R}^d} r_\alpha^2V(r)\,d^d\bm{r}=-\frac{1}{2d}\int_{\mathbb{R}^d} r^2V(r)\,d^d\bm{r}$

那么我们有

$\frac{1}{2}\int_{\mathbb{R}^d}\rho(\bm{r}) V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'=\gamma\rho^2(\bm{r})-\frac{\kappa}{2}\rho(\bm{r})\Delta\rho(\bm{r})+o(\nabla^4\rho)$

我们把一个卷积形式的项改写为了局域项。在最后积分时如果再对拉普拉斯项做一次分部积分，假设边界项为0，我们可以把自由能泛函写成

$\mathcal{F}[\rho]=\int_{\mathbb{R}^d}\left(V_{\mathrm{ext}}(\bm{r})\rho+\gamma\rho^2+\frac{\kappa}{2}(\nabla\rho)^2+k_BT\rho\ln\rho\right)\,d^d\bm{r}$

这正是标准Cahn-Hilliard自由能的形式

$\mathcal{F}[\rho]=\int_{\mathbb{R}^d}\left(g(\rho)+\frac{\kappa}{2}(\nabla\rho)^2\right)\,d^d\bm{r}$

因为所有的量都是局域的，所以自由能密度也成了 $\rho$ 的函数

$f(\rho)=V_{\mathrm{ext}}(\bm{r})\rho+\gamma\rho^2+\frac{\kappa}{2}(\nabla\rho)^2+k_BT\rho\ln\rho$

十一、相变和相平衡条件

我们可以对无外场时的动力学方程

$\frac{d\rho}{dt}=\nabla\cdot\frac{1}{\mu}\int d\bm{r}'\,\nabla_{\bm{r}} V(\bm{r}-\bm{r}')\rho(\bm{r'})\rho(\bm{r})+\frac{k_BT}{\mu}\nabla^2\rho$

进行线性稳定性分析，令 $\rho=\rho_0+\int_{\mathbb{R}^d}\delta\rho(\bm{k}) e^{i\bm{r}\cdot\bm{k}}\,d^d\bm{k}$ 保留到 $\delta\rho$ 的线性项，来判断粒子的均匀分布失稳的条件。相应的分岔条件给出的粒子密度 $\rho_0$ ，称为spinodal密度。在spinodal密度给出的区间内，粒子的均匀状态是不稳定的，也就是说此时系统内不存在单一相，系统必然会相分离——气体会凝结成液滴。简单计算可以得到失稳条件是

$k_BT+\rho_0V(k)<0$

然而这个条件只是线性失稳条件。完全相分离时，两相的密度我们称为binodal密度。非线性的作用使得spinodal和binodal并不重合。观察相图，往往能看到参数空间中的一个区域在binodal曲线的内部，而在spinodal曲线的外部。在这个区域内，系统的均匀状态是线性稳定的，但相分离的状态也是稳定的。大的扰动有可能使均匀的状态失稳，而造成相分离。这正类似于过冷水结冰或者过饱和水蒸汽凝结。

为了计算相分离时的两相密度，我们来寻找相分离条件。前一节实际上已经给出了一个条件：平衡态时化学势

$\nu=\frac{\delta\mathcal{F}[\rho]}{\delta\rho}=\int_{\mathbb{R}^d} V(\bm{r}-\bm{r}')\rho(\bm{r}')\,d^d\bm{r}'+k_BT\ln\rho(\bm{r})+k_BT$
是常数，这和热力学保持一致。但这一个方程还无法解出两相的密度。对于另一个条件，从热力学中我们可以猜出来是压强也是常数。为了理解这一点，我们我们要尝试找到一个应力张量 $\bm{\sigma}$ ，使得

$\rho\nabla\frac{\delta\mathcal{F}}{\delta\rho}=-\nabla\cdot\bm{\sigma}=-\sum_\beta\partial_\beta\sigma_{\alpha\beta}$

那么第二个相平衡条件即是 $\nabla\cdot\bm{\sigma}=0$ ，或者说应力张量为常数。之所以说它是应力，是因为如果给系统放一个外场，系统的粒子流可以写为

$\bm{J}=\frac{1}{\mu}(\nabla\cdot\bm{\sigma}-\rho\nabla V_{\mathrm{ext}})$

如果将外场受力在任意一个区域V内做体积分，那么受力转换为一个面积分

$-\int_V \rho\nabla V_{\mathrm{ext}}\,dV=-\int_V\nabla\cdot\bm{\sigma}\,dV=-\oint_{\partial V}\bm{\sigma}\cdot d\bm{S}$

于是我们看到， $\bm{\sigma}\cdot d\bm{S}$ 的确是粒子在 $d\bm{S}$ 面上提供的压强，与外力平衡，所以这个相平衡条件的确是压强相等，或者说力学平衡。这两个条件是独立的，因为一个相当于对 $\nabla(\delta\mathcal{F}/\delta\rho)$ 积分，另一个相当于对 $\rho\nabla(\delta\mathcal{F}/\delta\rho)$ 积分。

如果粒子相互作用的势是球对称的，经过一个较为复杂的计算（可以参考Exact formulation and coarse grained theory），我们可以得到

$\sigma_{\alpha\beta}=-k_BT\rho\delta_{\alpha\beta}+\frac{1}{2}\int d^d\bm{r'}\,\frac{r_\alpha' r_\beta'}{r'}\frac{dV(r')}{dr'}\int_0^1d\lambda\,\rho(\bm{r}-\lambda\bm{r}')\rho(\bm{r}+(1-\lambda)\bm{r}')$

第一项是熟悉的理想气体压强，它是各项同性的。而第二项则是相互作用带来的修正，容易理解它在均匀相内也是各项同性的。但是，在两相交界处，它就不再是各项同性的了。这时我们可以找到表面张力——也就是在气液界面处法向与切向的压强差。如果界面垂直于z轴，系统有xy平面内的平移对称性，考虑积分的两端在两相内部深处，那么我们可以定义表面张力为

$\begin{align} \gamma=&\int_{z_g}^{z_l}((\sigma_{xx}+\sigma_{yy})/2-\sigma_{zz})\,dz \\ =&\frac{1}{2S}\int d^d\bm{r}d^d\bm{r}'\, \frac{(x'^2+y'^2)/2-z'^2}{r'}\frac{dV(r')}{dr'}\int_0^1d\lambda\,\rho(\bm{r}-\lambda\bm{r}')\rho(\bm{r}+(1-\lambda)\bm{r}') \end{align}$

其中S是界面的面积。如果界面处弯曲，考虑一个球形的气泡，假设气泡半径为R，那么我们在极坐标系内沿径向做积分

$\int_{r_g}^{r_l}\nabla\cdot\bm{\sigma}\cdot d\bm{r}=\int_{r_g}^{r_l}\left(\partial_r\sigma_{rr}+\frac{1}{r}(\sigma_{rr}-\sigma_{\theta\theta})\right)\,dr=0$

考虑到积分几乎只在气泡界面 $r_f$ 处有贡献，有

$\sigma_{rr}(r_l)-\sigma_{rr}(r_g)=-\int_{r_g}^{r_l}\frac{1}{r}(\sigma_{rr}-\sigma_{\theta\theta})\,dr\approx\frac{\gamma}{r_f}$

即气泡内外侧的压强差与气泡大小成反比，与表面张力成正比，这正是拉普拉斯压强。至此我们发现这个系统的热力学和传统的热力学定律并没有什么区别。

在结束整个平衡态统计的内容之前，最后再提一下，对于类似Cahn-Hilliard的自由能泛函（比如朗道-金兹堡理论），一般地可以写成

$\mathcal{F}[\rho]=\int f(\rho,\nabla\rho)d^d\bm{r}$

这类自由能泛函没有卷积项，是局域的，那么应力张量可以简单地得到

$\sigma_{\alpha\beta}=-\left(\rho\frac{\delta\mathcal{F}}{\delta\rho}-f\right)\delta_{\alpha\beta}-\frac{\partial f}{\partial(\partial_\alpha\rho)}\partial_\beta\rho$

我们依然有压强项和表面张力项，与之前的分析并没有什么不同。

十二、化学反应系统

最后再用一点点篇幅简单介绍一下这一套方法在化学反应系统中的应用。我们考虑充分混合均匀的化学反应系统，从而不考虑空间坐标。这时我们的相空间变为了系统各种分子的分子数构成的空间 $\bm{x}=(x_1,\cdots,x_N)$ ，其中N是系统中分子的种类数。我们考虑有M个可能的反应，其中第j个反应发生的概率为指数分布，且发生率为 $a_j(\bm{x})$ ，发生这个反应后系统的状态从 $\bm{x}$ 变为 $\bm{x}+\bm{\nu}_j$ 。（比如 $2H_2+O_2=2H_2O$ ， $\bm{\nu}=(-2,-1,2)$ ），那么系统的状态作为随机过程，依然是马尔可夫的。显然这个系统的主方程是

$\frac{\partial p(\bm{x};t)}{\partial t}=\sum_{j=1}^{M}a_j(\bm{x}-\bm{\nu}_j)p(\bm{x}-\bm{\nu}_j;t)-\sum_{j=1}^{M}a_j(\bm{x})p(\bm{x};t)$

对于二元反应，比如， $x_1+x_2\to x_3$ ，考虑到反应是分子间短程相互作用——即碰撞的结果，有确定能量的分子在碰撞中应该有确定的概率反应，其概率由量子力学决定。那么反应发生的概率应该正比于粒子发生碰撞的概率，固定体系体积则正比于分子数，即 $a(x_1,x_2)=kx_1x_2$ 。一般来说k和系统体积有关系。如果我们计算系统的平均分子数

$\langle\bm{x}\rangle=\sum_\bm{x} \bm{x}p(\bm{x};t)$

那么它满足

$\frac{d\langle\bm{x}\rangle}{dt}=\sum_{j=1}^M \bm{\nu}_j\langle a_j(\bm{x})\rangle$

称为反应速率方程。但一般来说a(x)是非线性函数，这个方程不封闭，我们依然需要"平均场假设" $\langle a_j(\bm{x})\rangle=a_j(\langle \bm{x}\rangle)$ 。这意味着系统内分子数很大，两种分子的涨落之间的关联小到可以忽略。这时方程封闭，且如果对于d元反应有 $k\propto V^{-d}$ ，我们又得到了质量作用定律。

注意到，当粒子数极多的时候，单个反应对系统状态的改变很小，作泰勒展开保留至二阶项，主方程可以近似为偏微分方程

$\frac{\partial p(\bm{x};t)}{\partial t}=-\sum_{j=1}^{M}\sum_{i=1}^N\nu_{ji}\frac{\partial a_j(\bm{x})p(\bm{x};t)}{\partial x_i}+\frac{1}{2}\sum_{j=1}^{M} \sum_{k,l=1}^N\nu_{jk}\nu_{jl} \frac{\partial^2 a_j(\bm{x})p(\bm{x};t)}{\partial x_k\partial x_l}$

考虑到 $\nu$ 一般情况下都是常数，这其实正是Fokker-Planck方程的形式。对照之前的结果，也就是说这时候系统可以用化学朗之万方程描述

$\frac{d\bm{x}}{dt}=\sum_{j=1}^M a_j(\bm{x})\bm{\nu}_j+\sum_{j=1}^{M}\sqrt{a_j(\bm{x})}\bm{\nu}_j\bm{\xi}_j$

其中 $\bm{\xi}_j$ 是delta关联的高斯白噪声。

最后简单提一下，从原始的主方程出发硬算，我们可以得到化学反应系统中的涨落耗散定理

$\frac{d\sigma_{kl}}{dt}=\langle A_k(\bm{x})(x_l-\langle x_l\rangle)+A_l(\bm{x})(x_k-\langle x_k\rangle)\rangle+\langle B_{kl}(\bm{x})\rangle$

其中 $\sigma_{kl}=\langle(x_k-\langle x_k\rangle)(x_l-\langle x_l\rangle)\rangle$ 是协方差矩阵，而

$A_k(\bm{x})=\sum_{j=1}^M \nu_{jk}a_j(\bm{x})$

$B_{kl}(\bm{x})=\sum_{j=1}^M \nu_{jk}\nu_{jl}a_j(\bm{x})$

正是在Fokker-Planck方程近似中出现的系数

$\frac{\partial p(\bm{x};t)}{\partial t}=-\sum_{i=1}^N\frac{\partial A_i(\bm{x})p(\bm{x};t)}{\partial x_i}+\frac{1}{2}\sum_{k,l=1}^N \frac{\partial^2 B_{kl}(\bm{x})p(\bm{x};t)}{\partial x_k\partial x_l}$

也就是说，我们再一次将平衡系统的涨落项与耗散项联系了起来。

来源：知乎 www.zhihu.com
作者：赵永峰

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

经典网文分享

软物质的统计力学（完）

目录

十、密度泛函理论

十一、相变和相平衡条件

十二、化学反应系统