温度涨落
想不到，教科书里的内容也有很大讨论的余地。随着分子动力学凝固形核模拟的发展，这件事必须再拿出来说一说了。本帖子是一篇文献的读书笔记：
Grey Sh. Boltachev and Jürn W. P. Schmelzer, J. Chem. Phys. 133, 134509 (2010); doi: 10.1063/1.3486557
传承了毛子一贯的厚重作风，从朗道一直讨论到21世纪。
基于爱因斯坦“涨落的热力学理论”的简单推导
整个系统的setup如图。整个系统是一个孤立系统。其性质用带下标tot的量表示。
假设系统的内能发生了改变，此改变不仅包含体积功。根据整个系统是一个孤立系统的假设，有：
\Delta U = W_{s}+P_0 \Delta V_0 - T_0 \Delta S_0\quad (1)
整个大系统的体积不变，所以有：
\begin{equation} \Delta V + \Delta V_0 = 0\quad (2) \end{equation}
熵变：
\Delta S + \Delta S_0 \ge 0\quad (3)
把这俩条件(2), (3)带入(1)，有：
W_s \ge \Delta U + P_0 \Delta V - T_0 \Delta S\quad (4)
则最小非体积功为：
{\rm{d}}W_{\rm{min}} = {\rm{d}}U + P_0 {\rm{d}}V - T_0 {\rm{d}}S = \Delta(U+P_0 V - T_0 S)\quad(5) （这里 P_0 和 T_0 可视为常数）
根据热力学基本定理： {\rm{d}}U = T{\rm{d}}S - P {\rm{d}}V
有：
{\rm{d}}W_{\rm{min}} = (T- T_0) {\rm{d}}S-( P-P_0) {\rm{d}}V \quad(6)
根据Landau卷四第二章$20的论述，可见这个最小非体积功，与平衡时的大系统的熵 S_{\rm{tot}}(U_{\rm{tot}}) 、系统与热源分别平衡，但互相未达到平衡时的熵之和 S_{\rm{tot}} 与 S_{\rm{tot}}(U_{\rm{tot}}) 之差有上图中的关系。极限情况下，就有：
\Delta S_{\rm{tot}}= -\frac{{\rm{d}}S}{{\rm{d}}U}{\rm{d}}W_{\rm{min}} = -\frac{W_{\rm{min}}}{T_0}\quad (7)
\Delta S_{\rm{tot}}= -\frac{\Delta (U +P_0 V - T_0 S)}{T_0} （这里 P_0 和 T_0 可视为常数）
根据Boltzmann熵定理，有：
\Delta S_{\rm{tot}} = k_{\rm{B}}\ln \frac{\Omega }{\Omega_0}
所以可以推出，当系统与环境分别达到热平衡，但两者有一定差别，这个差别非常小的时候，出现差别的概率：
{\rm{Probability}} \propto e^{-\frac{\Delta U + P_0\Delta V - T_0 \Delta S}{k_{\rm{B}}T_0}} \quad (8)
涨落很小的时候，可以展开 U = U(S,V) ，根据上一篇帖子和知友 @白白 的计算：
{\rm{Probability}}\propto e^{-\frac{(\Delta T \Delta S - \Delta P \Delta V)}{2k_{\rm{B}}T_0}}
注意这里跟上边俩帖子不同，考虑了系统与环境未达到平衡时的情况。
再选取 T 、 V 为自变量，可以得到：
{\rm{Probability}}\propto \exp\left\{ -\alpha (\Delta T/T_0)^2 -\beta(\Delta V/V_0)^2 \right\}
\alpha = \frac{C_{V}}{2k_{\rm{B}}},\quad \beta = -\frac{V^{2}_0}{2k_{\rm{B}}T_0}\left( \frac{\partial P}{\partial V} \right)_{T_0}
这样保证指数上无量纲，积分限也不带量纲，方便下边推导。
取热力学极限的积分
这里，如果积分限可以取为正负无穷，有：
\langle \Delta T\rangle = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\Delta T e^{-\alpha (\Delta T)^2-\beta (\Delta V)^2}{\rm{d}}(\Delta T/T_0){\rm{d}}(\Delta V/V_0)}{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}e^{-\alpha (\Delta T)^2-\beta( \Delta V)^2}{\rm{d}}(\Delta T/T_0){\rm{d}}(\Delta V/V_0)} = 0
而且，
\langle T \rangle = T_0,\quad \langle (\Delta T)^2\rangle = \frac{k_{\rm{B}}T^2_0}{C_V}
就是课本上的结论。Landau的书讲的是最仔细的了。
非热力学极限的积分
但故事并不这样简单，因为积分限实际上不能达到负无穷。积分下限应该是 \Delta T = -T_{0},\quad \Delta T/T_0 = -1 。
\langle \Delta T\rangle = \frac{\int_{-1}^{\infty} \int_{-\infty}^{\infty}\Delta T e^{-\alpha (\Delta T)^2-\beta (\Delta V)^2}{\rm{d}}(\Delta T/T_0){\rm{d}}(\Delta V/V_0)}{\int_{-1}^{\infty} \int_{-\infty}^{\infty}e^{-\alpha (\Delta T)^2-\beta( \Delta V)^2}{\rm{d}}(\Delta T/T_0){\rm{d}}(\Delta V/V_0)}
\langle T\rangle = T_0\left\{1+ \frac{k_{\rm{B}}}{C_{V}}\frac{\exp\left(-\frac{C_{V}}{2k_{\rm{B}}}\right)}{ \int_{-1}^{\infty} e^{-\alpha x^2}{\rm{d}}x } \right\}
这时候会发现，系统的平均温度要高于环境的。显然，在热力学极限下， \langle T \rangle = T_0 ；但是如果考虑一个小系统，就会出现这种情况。而且，如果 N\to \infty ，这个积分的被积函数会尖锐的集中在 \Delta T = 0 附近，所以对积分限的任何延拓都是可以的。但一旦体系是少原子的情况，根据上述“标准分析”就会出现这个矛盾。
解决矛盾
为了解决上述问题，重新回到最小功的公式（6），考虑恒容情况：
{\rm{d}}W_{\rm{min}} = (T-T_0){\rm{d}}S
考虑到：
C_{V}=T \left(\frac{\partial S}{\partial T}\right)_{V}
得到：
{\rm{d}}W_{\rm{min}} = C_{V}\frac{T-T_0}{T}{\rm{d}}T
假设热容不随温度变化，积分此等式：
W_{\rm{min}} = C_{V}\int_{T_0}^{T}(1-\frac{T_0}{T}){\rm{d}}T
W_{\rm{min}} = C_{V}\left[(T-T_0)-T_0\ln\frac{T}{T_0}\right] （如果热容是一般情况，需要通过分子动力学模拟得到热容，再做数值积分，没那么简单。）
根据Einstein涨落公式（整不明白就参考朗道第四卷吧）：
{\rm{Probability}} = A_1 e^{-\frac{W_{\rm{min}}}{k_{\rm{B}}T_0}}
根据这个公式，可以证明：
\langle \frac{1}{T}\rangle = \frac{1}{T_0}
\langle T\rangle =T_0(1+\frac{k_{\rm{B}}}{C_{V}})
非常神奇。
这里可以看到，如果一个系统足够小，其温度就跟热源不一样，这样岂不是可以设计一个第二类永动机了？这就是为什么C. Kittel不承认存在温度涨落的原因。因为只有取到热力学极限，才不会发生这种佯谬。这也是热物理学的本意：只讨论大量（宏观意义上）粒子的系统。只基于系综理论，讨论能量涨落，不讨论温度的涨落就没事了。
正反两方的意见在此：
- C. Kittel, Phys. Today 41 5 , 93 1988 .
- B. B. Mandelbrot, Phys. Today 42 1 , 71 1989 .
这个问题会使凝固形核的机理解释不清楚：
- J. Wedekind, D. Reguera, and R. Strey, J. Chem. Phys. 127, 064501 2007 .
- J. Wedekind, D. Reguera, and R. Strey, in The Temperature of Nucleating Droplets, Nucleation and Atmospheric Aerosols, edited by C. D. O’Dowd and P. E. Wagner Springer, New York, 2007 , pp. 102–106.
Schmeltzer提出解决之道：
考虑单一原子的理想气体，由理想气体的内能公式：
U = \sum_{i=1}^{N}\frac{m v^2_i}{2}\quad (9)
以及能均分定理：
U = \frac{3Nk_{\rm{B}}T}{2}
得：
T=\frac{m}{3k_{\rm{B}}N}\sum_{i=1}^{N}v^{2}_i\quad(10)
公式（10）正是做MD模拟用到的计算系统温度的公式。做过NVT系综的MD模拟，就会深入理解“温度的涨落”和“温度的计算”有多么重要。
再根据理想气体的Maxwell-Boltzmann分布律：
f(v) = \prod_{i=1}^{N}f(v_i) = A\prod_{i=1}^{N}v^2_{i}e^{-\frac{m\beta_0}{2} v^2_i}
这里 A 是归一化系数。
利用(9)，把上边的M-B分布改为能量分布（利用N维球体积公式和表面积公式，把速度都积分掉，最后再做对能量的积分把系数定出来）：
f(U) = \frac{1}{\Gamma(3N/2)}\frac{1}{U}\left(\frac{U}{k_{\rm{B}}T_0}\right)^{3N/2} e^{- \beta_0 U}
注意这里温度是热源的。
这时候，由于内能有确定的瞬时值，所以可以定义瞬时温度的分布表达式（注意归一化）：
f(T) = \frac{1}{\Gamma(3N/2)}\frac{1}{T}\left(\frac{3NT}{2 T_0}\right)^{3N/2}e^{-3NT/(2 T_0)} \quad (11)
用这个瞬时温度的分布表达式，可以证明其平均值就是热源的温度 T_0 ：
\langle T \rangle = \int_0^{\infty}Tf(T){\rm{d}}T = T_0
此时对理想气体，其平方平均：
\begin{aligned} \langle T^2 \rangle & = \int_0^{\infty} T^2 f(T) {\rm{d}}T = (1+\frac{2}{3N})T^2_0\\ \end{aligned}
其涨落仍然接近 \frac{1}{\sqrt{N}}
这个瞬时温度的分布表达式，跟温度的准静态理论有很大区别。根据朗道的书里的表达式，代入理想气体的相关信息，是：
f(T) = \sqrt{\frac{3N}{4\pi T^2_0}}e^{-\frac{3N}{4}\left(\frac{T-T_0}{T_0}\right)^2},\quad(C_V = \frac{3N}{2})
实际上这就是Gamma-分布和Gauss分布的区别。根据概率论，在 N 很大的时候两者相同，但小 N 时差别很大。
文章作者提出，对于真实流体，可以用
C_{V}=\frac{3Nk_{\rm{B}}}{2}
代换掉粒子数 N ：
f(T)_{\rm{MD}} = \frac{1}{\Gamma(C_{V}/(k_{\rm{B}}))}\frac{1}{T}\left( \frac{C_{V}T}{k_{\rm{B}}T_{0}} \right)^{C_{V}/k_{\rm{B}}}\exp\left(-\frac{C_{V}T}{k_{\rm{B}}T_{0}}\right)
后记
而利用MD计算 C_{V} ，可根据定义直接计算，或者拟合插值得到，又是另一个坑了。
由公式（35）推导公式（40）：
\begin{aligned} \frac{1}{T}=\frac{1}{T_0} - \frac{1}{C_{V}T}\left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right)\end{aligned}\quad (35)
\begin{aligned}1 = \frac{T}{T_0} - \frac{T}{C_{V}T_0}\left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right)\end{aligned}
\begin{aligned} T_0 = T - \frac{T}{C_{V}} \left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right) \end{aligned}
\begin{aligned} T = T_0 + \frac{T}{C_{V}}\left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right) \end{aligned} \quad ({\rm{W}})
把公式(W)代入（37）：
\begin{aligned} \langle T \rangle =\langle \left[ (T_0 + \frac{T}{C_{V}}\left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right)\right]\rangle \end{aligned}
这里尖括号代表求热力学平均值：
\langle g(T) \rangle = \frac{\int_{0}^{\infty}g(T)e^{-W_m/(k_{\rm{B}}T_0)}{\rm{d}}T}{\int_{0}^{\infty}e^{-W_m/(k_{\rm{B}}T_0)}{\rm{d}}T}
右边第二式可以分部积分：
\begin{aligned} & \langle \frac{T}{C_{V}}\left(\frac{{\rm{d}}W_{m}}{{\rm{d}}T}\right)\rangle \\ = &\frac{1}{\int_{0}^{\infty}e^{-W_m/(k_{\rm{B}}T_0)}{\rm{d}}T} \int_{0}^{\infty} \frac{T}{C_{V}} e^{-W_{m}/(k_{\rm{B}}T_0)}{\rm{d}}W_{m}\\ = & \frac{k_{\rm{B}}T_0}{C_{V}\int_{0}^{\infty}e^{-W_m/(k_{\rm{B}}T_0)}{\rm{d}}T}\int_{0}^{\infty} T e^{-y} {\rm{d}}y \\ = & \frac{k_{\rm{B}}T_0}{C_{V}\int_{0}^{\infty}e^{-W_m/(k_{\rm{B}}T_0)}{\rm{d}}T}\int_{0}^{\infty} e^{-W_{m}/(k_{\rm{B}}T_0)}{\rm{d}}T \\ = & \frac{k_{\rm{B}}T_0}{C_{V}} \end{aligned}
所以有：
\langle T\rangle = T_0 \left(1+ \frac{k_{\rm{B}}}{C_{V}}\right)
最后推荐一下朗道的统计物理学。这本书第一章非常提纲挈领，需要把整个一门“热统”学完了才能看懂。所以这本书不建议第一次学习，纯属研究生级别的课本。但非常juicy，时常能读出新想法。