由最小作用量原理可知，能量项体现了物理系统的本质特征，这一观点同样适用于平衡态的统计系统。在统计物理中，平衡态的配分函数正是建立在能量项基础之上的。此处涉及能量零点的选取问题。由于统计系统的温度反映的是粒子热运动动能的平均值，因此需排除介质整体平动所带来的影响。换言之，配分函数的变量应为共动系中的热运动能量项。于是，介质的能量分布函数可写为$dP=f(K/kT)dK$，其中$K$表示热运动动能，是随机变量；$T$为系统温度；$k=1.38\times 10^{-23}$ J/K是玻尔兹曼常数。该分布函数满足以下矩方程：

 $$ \int _0^\infty f \left(\frac K{kT}\right) d K=1,\quad \int _0^\infty K d {P}=\frac 3 2 kT,\quad \int _0^\infty K^2 d P=\frac 3 {2\sigma} (kT)^2, \qquad  (1) $$

其中，第二个公式定义了介质的温度$T$，$\sigma$是刻画粒子动能分布的常数。在理想气体的麦克斯韦分布下，$\sigma=2/5$。

 在热力学实验中，常用气缸活塞装置研究气体的物态方程。然而，该装置仅为近似模型：由于活塞运动引发气体的宏观流动，破坏了平衡态，只有在无穷小速度下才能逼近理想过程。但在宇宙尺度上，这一问题可以克服——时空本身便是理想的气缸-活塞装置，因为这是一个共动参考系。在大尺度下，时空具有均匀性和各向同性，符合FLRW度规： 

$$ ds^2=d\tau^2 -a(\tau)^2\left(dr ^2+S (r)^2(d \theta^2+\sin^2\theta d\varphi^2)\right),$$

其中$S(r) =\sin r, r, \sinh r$，但下面的讨论对$S(r)$具体形式并无依赖关系。为简单起见取$c=1$作为速度单位。体积元与$a^3$成正比，$a$的变化即代表“活塞”运动。由于该度规具备高度对称性，对于所有物质，物质拉格朗日量${\cal L}_m$产生的能量-动量张量具有传统形式：

$$ T^\mu_\nu={\rm diag}(\rho,-P,-P,-P).  $$

根据能量守恒或热力学第一定律，有

$$ T^{\mu\nu}_{;\nu}=0 \Leftrightarrow d (\rho V) + P d V =0 \Leftrightarrow d({\rho a^3})=-P da^3. \qquad (2)  $$

对于FLRW度规，唯一独立的动力学约束即为Friedmann方程：

$$ \dot a^2 =\frac 1 3 \Lambda a^2-K + \frac \kappa 3\rho a^2.  $$

给定初始条件与物态方程$P=P(\rho)$，上述方程构成封闭系统，存在唯一解[1]。

我们首先分析稀薄理想气体的情况，涵盖单原子气体和光子；其他情形可按类似逻辑处理。此时每个粒子的动量满足：

$$ p_n(\tau) a(\tau)=p_n(\tau_0)a(\tau_0)=C_n,\qquad p= \frac {m_n v_n}{\sqrt {1-v_n^2}}, \qquad (3) $$

其中$m_n$为第$n$个粒子的固有质量。对无质量光子可证，其波长$\lambda(\tau)\propto a(\tau)$，从而动量$p$亦满足(3)。

由此建立理想气体的状态函数$(\rho, P)$与尺度因子$a$之间的关系。由(3)得 $p^2_n={C_n}^2{a^{-2}}$，其中$C_n$由初始时刻$\tau_0$的数据确定。因此，对所有粒子，其动量平方平均值与参数$a$的关系为 $\bar p^2= {C_0}a^{-2}$，$C_0$为常数。

另一方面，相对论下动量$p$与动能$K$关系为：

$$ p^2= K( K+2m),\qquad K=\frac m{\sqrt{1-v^2}}-m, \qquad (4) $$

由(4)结合(1)得：

$$ \bar p^2=\frac{1}{N}\sum_{n=1}^N \int^\infty_0 p^2_n f( K_n)d K_n = \frac {3}{2\sigma} kT(kT+{2\sigma}\bar m), \qquad (5)  $$

其中$\bar m=\frac 1 N \sum_n m_n$为平均静质量。比较(5)与$\bar p^2= {C_0}a^{-2}$，可得温度$T$与$a$的关系：

$$ kT=\frac{\sigma\bar m} a (\sqrt{a^2+{b^2}}-a),\qquad a=\frac {\sigma b\bar m}{\sqrt{kT(kT+2\sigma\bar m)}}, \qquad  (6) $$

其中常数$b>0$由初始条件确定。

在后续统计计算中，以$a$为自变量。由(6)得统计平均值：

$$ \rho  \equiv  \sum_n \frac {m_n}{\sqrt{1-v^2_n}}  \delta^3(\vec x -\vec X_n) = \frac 1 { V} \int_{ V} \sum_{X_n\in V}\frac {m_n } {\sqrt{1-v_n^2}}  \delta^3(\vec x-\vec X_n)\sqrt{\bar g} d^3x $$

$$  =  \frac 1 {a^3  \Omega}\sum_{X_n\in  \Omega}\int\frac {m_n } {\sqrt{1-v_n^2}}d P =\frac 1 {a^3  \Omega} \sum_{X_n\in  \Omega} \int (m_n+ K_n)d P $$

$$  = \frac 1 {a^3  \Omega}\sum_{X_n\in  \Omega} (\bar m+\frac 3 2 kT)=\frac {\varrho}{a^3} \left( 1+\frac {3\sigma} {2a} {(\sqrt{a^2+b^2}-a)}\right),\qquad (7) $$

其中$ V=a^3 \Omega$, $\varrho =\frac 1  \Omega \sum_{X_n\in  \Omega}m_n$是静质量的角密度，这是一个与$a$无关的常数。将(7)代入能量守恒定律(2)可得

$$ P \equiv \frac 1 3 \sum_n \frac {m_n v^2_n}{\sqrt{1-v^2_n}}  \delta^3(\vec x -\vec X_n)=\frac {\sigma\varrho b^2}{2a^4\sqrt{a^2+b^2}}=\frac {NkT} { V} \left(1-\frac{ kT}{2(\sigma\bar m+kT)}\right). $$

当$a\to 0$时，$\rho\sim P\to C a^{-4}$，说明理想气体无法提供负压。

上述推导将FLRW度规视为驱动理想气体的活塞。微观角度看，粒子与光子在平均引力作用下沿测地线运动，粒子间的碰撞近似为瞬时事件。由此，所有热力学量均可经动力学与统计学精确求解。假设$N\gg 1$个粒子分布于体积微元$\Delta V$内，此即“爱因斯坦电梯”模型，但此电梯非刚体，其体积在引力作用下随时间变化。由于热力学方程（如状态方程、能量守恒）为局部统计规律，与宏观几何形状无关，因此采用对称的FLRW时空推导状态方程具备简洁明晰的优势。根据等效原理，这些结果在共动系中普遍适用。进一步分析可知，粒子的弹性碰撞对结论无影响，所以上述推论具有普适性[2]。由(6)得： 

$$ a = \frac {\sigma b}{\sqrt{J(J+2\sigma)}},\qquad J \equiv \frac{kT}{\bar m c^2}=\frac\sigma a (\sqrt{a^2+b^2}-a), $$ 

其中$J$为无量纲温度。

综合以上结果得到如下结论：

定理：对于相对论性理想气体，物态方程为： 

$$ P=\frac {\bar\rho kT} {\bar m c^2} \left(1-\frac{ kT}{2(\sigma\bar m c^2+kT)}\right),\quad \bar m \equiv \frac 1 N \sum_{n=1}^N m_n,\quad \bar\rho\equiv \frac {N\bar m} V\propto \frac 1 {a^3}. $$ 对于绝热过程，状态函数满足参数方程： $$ \bar\rho = \rho_0 [J(J+2\sigma)]^{\frac 3 2},\quad \rho = \bar \rho (1+\frac 3 2 J),\quad P = \bar\rho \frac {J(J+2\sigma)} {2(J+\sigma)},  $$ 其中$J$为参数，$\rho_0= \varrho{(\sigma b)^{-3}}$是带密度量纲的常数。

上述推导满足洛伦兹不变性，且考虑了引力场的作用，结论与相对论是相容的，因此这个物态方程可以避免黑洞产生。在低温极限下，定理给出单原子气体的绝热状态方程：

 $$ P\dot = \rho J \dot = P_0 \rho^{\frac 5 3},\qquad (J\ll 1,~{\rm or}~kT\ll \bar m),  $$ 

这与实验热力学定律相符。当$J\to \infty$或$\bar m\to 0$时，有$P\to \frac 1 3 \rho$，并得Stefan–Boltzmann定律$\rho \propto T^4$，因此该模型亦适用于光子或极端相对论性粒子。在广义相对论中，所有过程自然发生，常数$\rho_0$仅由气体性质决定，与具体过程无关。

以上推导基于点粒子的动力学模型：

$$  {\mathcal L} ={\mathcal L}_g+{\mathcal L}_p,\quad {\mathcal L}_g =\frac1 {16\pi G}(R-2\Lambda),  $$

其中${\mathcal L}_p$表示无相互作用点粒子的拉格朗日量：

$$  {\mathcal L}_p = -\sum_n m_n \sqrt{1-v_n^2}\delta^3(\vec x-\vec X_n),\quad v^2_n=g^{00}\bar g{kl}\frac {dX_n^k}{dt}\frac {d X_n^l}{dt}. \qquad (8)  $$

该点粒子模型涉及 $\delta^3(\vec x - \vec X_n)$，这一表达只有在自然坐标系下才具有明确定义：

$$  d{\bf x}^2 = g_{00} dt^2-\bar g_{kl}dx^k dx^l,\quad \sqrt{g}=\sqrt{g_{00}}\sqrt{\bar g},\quad {\bar g}={\det(\bar g_{kl})}, \qquad (9)  $$

其中，$d\tau = \sqrt{g_{00}} dt$ 类似牛顿意义下的绝对时间，它不同于粒子固有的本征时间 $ds_k = \sqrt{1 - v_k^2} d\tau$。自然坐标系作为一种客观存在，使得式 (8) 在广义坐标变换下具有一般性，从而保证了上述物态方程在相对论意义下的正确性；事实上，(8) 是洛伦兹不变量。统计平均是相对于气体微元的共动坐标系进行的，因此相关物理量如 $\rho, P, T$ 都应理解为固有量。在洛伦兹变换下，$T$ 是标量，而 $\rho$ 与 $P$ 的变换性质需额外处理，下面我们具体讨论这一问题。

若理想气体处于某一标量作用势 $\Phi$ 之中，其微观拉格朗日量为：

$$  {\mathcal L}_f=\sum_n (-m_n\sqrt{1-v_n^2}-\Phi)\delta(\vec x-\vec X_n), \qquad (10) $$

其中 $\Phi$ 表示相互作用势，虽然实际情况可能非常复杂，但通常处理较为直接。接下来我们探讨如何将式 (10) 简化为合理的宏观形式。设流体相对于自然坐标系 (9) 有一平均速度场 $\vec V$，在流体的共动坐标系 $\bar x^\mu$ 中，存在如下局部洛伦兹变换：

$$  \sqrt{1-v_n^2}=\sqrt{1-\bar v_n^2} \frac {d\bar t}{d\tau}=\sqrt{1-\bar v_n^2} \frac {\sqrt{1-\vec V^2}}{1-|\vec V|\bar v_n\cos\theta}, $$

其中 $\vec V^2 = \bar g_{kl} V^k V^l$，且 $\vec V \cdot \vec v_n = |\vec V|\bar v_n \cos \theta$。对于式 (7)所定义的质能密度，在统计意义下有：

$$\rho_p = \sum_n \frac {m_n}{\sqrt{1-v_n^2}}\delta(\vec x -\vec X_n)= \frac 1 {\sqrt{1-\vec V^2}}\sum_n \frac{m_n(1-|\vec V|\bar v_n\cos\theta)} {\sqrt{1-\bar v_n^2}}\delta(\vec x -\vec X_n) $$

$$\ =\frac 1 {\sqrt{1-\vec V^2}}\sum_n {m_n} \left(\frac 1 {4\pi}\int_\Omega \frac {(1-|\vec V|\bar v_n\cos\theta)\sin\theta d\theta d\varphi} {\sqrt{1-\bar v_n^2}} \right)\delta(\vec x -\vec X_n) $$

$$\ =\frac 1 {\sqrt{1-\vec V^2}} \sum_n \frac {m_n}{\sqrt{1-\bar v_n^2}} \delta(\vec x -\vec X_n) = \frac \rho {\sqrt{1-\vec V^2}}, \qquad \qquad $$

其中 $\rho$ 为共动坐标系中的质能密度，包含热能，$\Omega$ 表示单位球面。显然，该公式正确描述了以速度 $\vec V$ 移动的天体在自然坐标系中的质能密度。

基于 (8)，我们对 ${\mathcal L}_p$ 的平均可表示为：

$$ -{\mathcal L}_p = \sum_n m_n \sqrt{1-v_n^2}\delta(\vec x -\vec X_n)= \sqrt{1-\vec V^2}\sum_n m_n \frac {\sqrt{1-\bar v_n^2}}{1-|\vec V|\bar v_n\cos\theta}\delta(\vec x -\vec X_n) $$

$$\ = \sqrt{1-\vec V^2}\sum_n m_n \left(\frac 1 {4\pi}\int_\Omega \frac {\sqrt{1-\bar v_n^2}\sin\theta d\theta d\varphi}{1-|\vec V|\bar v_n\cos\theta} \right)\delta(\vec x -\vec X_n) $$

$$\ = \sqrt{1-\vec V^2}\sum_n {m_n} \frac {\sqrt{1-\bar v_n^2}}{|\vec V|\bar v_n} \ln\sqrt{\frac {1+|\vec V|\bar v_n}{1-|\vec V|\bar v_n}} \delta(\vec x -\vec X_n).$$

$$ \  \to  \sqrt{1-\vec V^2}\sum_n {m_n} {\sqrt{1-\bar v_n^2}} \delta(\vec x -\vec X_n)= \sqrt{1-\vec V^2}\rho-s P, $$

其中 $s$ 为待定常数，因为拉格朗日量是不同能量项的线性组合。

为获得正确的流体动力学方程，可证明 $s=1$，因此式 (10) 对应的宏观表达为： 

$$  {\mathcal L}_f =-\rho \sqrt{1-\vec V^2}+ P-\Phi. \qquad (11)  $$

在 (11) 中，自变量应为粒子位置 $\vec X$ 而非速度 $\vec V$，这也是流体拉格朗日量难以直接写出的主要原因之一。虽然该形式为现实流体提供了一个合理模型，但它无法导出连续性方程。此外，式 (11) 中的 $\rho$ 与 $P$ 含义不同于广义相对论中能量-动量张量的标准形式：

$$  T^{\mu\nu}\to (\rho+P) U^\mu U^\nu -P g^{\mu\nu},\qquad (12)  $$

对应(12)的流体方程 $T^{\mu\nu}_{\ ;\nu}=0$ 没有合适的拉格朗日量。根据 (11)，流体的能量-动量张量计算如下：

$$ T_f^{\mu\nu} =-2\frac{\delta ({\mathcal L}_f\sqrt{g})}{\sqrt{g}\delta g_{\mu\nu}}=\left\{\begin{array}{l}\rho\sqrt{1-V^2}{U^0 U^0} +(\Phi- P)g^{00},\\ \rho\sqrt{1-V^2} ({U^k U^l}+g^{kl}) +(\Phi- P)g^{kl}.\end{array}\right. \qquad (13)  $$

式 (13) 和 (12) 哪一种更为合理目前尚不明确，有待专门研究。在某些情况下，近似处理可能导致某些对称性破缺，例如经典近似和统计平均等，只有基本的旋量场和矢量场才具有完整的对称性。

在费米子量子统计中，通常使用 $e_n = \sqrt{p_n^2 + m_n^2}$ 作为随机变量。但值得注意的是，这里 $m_n$ 是常数，并非变量，只有动能 $K_n$ 才能定义温度 $T$，而 $e_n$ 不能。此外，直接以 $e_n$ 为基础进行能量分布估计，在推导白矮星的压强与质量上限时，可能会忽略关键的动力学效应，其结论是值得认真检查的。

参考文献

[1] Gu, Y.-Q. (2024) Why the Big Bang Never Happened[J]. Mathematics and Systems Science, 2 (2), 2630.(中文版)

[2] Functions of State for Spinor Gas in General Relativity, arXiv:0711.1243