玻色统计和费米统计

• 热力学统计量表达式
• 弱简并理想玻色气体和费米气体
• 玻色-爱因斯坦凝聚(BEC)
• 热辐射和光子气体
  • 瑞利-金斯公式
  • 维恩公式
  • 普朗克公式
• 金属中的自由电子气体

满足定域或非简并条件的粒子系统可以采用玻尔兹曼统计，而不满足非简并条件的系统，需要用玻色-爱因斯坦统计或费米-狄拉克统计处理。

不满足非简并条件的气体叫做简并气体。不同于理想气体，简并气体不满足定域条件。尽管气体之间的距离更大，但是气体粒子本身的波包尺度也更大。同样的道理，固体虽然原子之间间距小，但是波包更小，是定域系统。

粒子的简并情况也有所区别，具体可根据逸度$z=e^{-\alpha}$进行分类：

理想气体（非简并）	理想玻色/费米气体（弱简并）	玻色-爱因斯坦凝聚	光子气体	自由电子气（强简并）
$e^{-\alpha}\ll1$	$e^{-\alpha}\ll1$但不可忽略	$e^{-\alpha}\to 1^{-}$	$e^{-\alpha}= 1$	$e^{-\alpha}\gg 1$

其中，理想气体因为忽略相互作用，所以化学势$\mu$接近0，但我们在下面会讨论如果不忽略其影响会带来什么变化。自由电子气由于其费米能（化学势）远大于$kT$，所以在统计上表现出强简并特性。

我们多次谈到逸度，所以逸度是什么呢？由$e^{-\alpha}=e^{\mu/kT}$可得，逸度与化学势正比——当逸度很小的时候，意味着化学势远低于0，此时粒子倾向于逃逸出系统，这就是逸度的来源。当逸度很大的时候，化学势；远高于0，则会出现强简并的情况。

热力学统计量表达式

考虑玻色系统，对于巨正则系综，定义巨配分函数：

$$\Xi = \prod_l \Xi_l=\prod_l (1-e^{-\alpha-\beta\epsilon_l})^{-\omega_l}$$

在系综理论中，我们会推导玻色统计和费米统计的巨配分函数。

系综的平均粒子数：

$$\langle N \rangle = -\frac{\partial }{\partial \alpha} \ln \Xi=\sum_l \omega_l\frac{\partial }{\partial \alpha}\ln [1-e^{-\alpha-\beta\epsilon_l}] = \sum_l \frac{\omega_l}{e^{\alpha+\beta\epsilon_l}-1}$$

系综和平均能量：

$$\langle E \rangle = -\frac{\partial }{\partial \beta} \ln \Xi=\sum_l \omega_l\frac{\partial }{\partial \beta}\ln [1-e^{-\alpha-\beta\epsilon_l}] = \sum_l \frac{\omega_l\epsilon_l}{e^{\alpha+\beta\epsilon_l}-1}$$

外界的广义作用力

$$\langle Y \rangle = -\frac1\beta\frac{\partial }{\partial y} \ln \Xi=\sum_l \omega_l\frac{\partial }{\partial y}\ln [1-e^{-\alpha-\beta\epsilon_l}] = \sum_l \frac{\omega_l}{e^{\alpha+\beta\epsilon_l}-1}\frac{\partial \epsilon_l}{\partial y}$$

熵：

$$S =k(\ln \Xi + \alpha \langle N \rangle + \beta \langle E \rangle)$$

巨热力学势：

$$J = -kT\ln \Xi$$

对于费米系统，巨配分函数：

$$\Xi = \prod_l \Xi_l=\prod_l (1+e^{-\alpha-\beta\epsilon_l})^{\omega_l}$$

将玻尔兹曼统计和玻色统计和费米统计对比:

	玻尔兹曼统计	玻色统计	费米统计
配分函数	$Z_1 = \sum_l \omega_le^{-\beta\epsilon_l}$	$\Xi = \prod_l (1-e^{-\alpha-\beta\epsilon_l})^{-\omega_l}$	$\Xi = \prod_l (1+e^{-\alpha-\beta\epsilon_l})^{\omega_l}$
平均粒子数	$N = e^{-\alpha}Z_1$	$\langle N \rangle =-\frac{\partial }{\partial \alpha} \ln \Xi$	$\langle N \rangle =-\frac{\partial }{\partial \alpha} \ln \Xi$
平均内能	$U =N(-\frac{\partial}{\partial\beta})\ln Z_1$	$\langle U\rangle =-\frac{\partial }{\partial \beta} \ln \Xi$	$\langle U \rangle = -\frac{\partial }{\partial \beta} \ln \Xi$
外力	$Y -\frac{N}{\beta}\frac{\partial }{\partial y}\ln Z_1$	$\langle Y \rangle =-\frac1\beta\frac{\partial }{\partial y} \ln \Xi$	$\langle Y \rangle =-\frac1\beta\frac{\partial }{\partial y} \ln \Xi$
熵	$S=Nk(\ln Z_1+\beta U)$	$S =k(\ln \Xi + \alpha \langle N \rangle + \beta \langle E \rangle)$	$S =k(\ln \Xi + \alpha \langle N \rangle + \beta \langle E \rangle)$

其实，玻色和费米统计的配分函数比玻尔兹曼“高一阶”（从求和变成连乘可以看出端倪）。对其求对数：

$$\ln \Xi = \sum_l \ln (1\mp e^{-\alpha-\beta\epsilon_l})^{\mp\omega_l} = \sum_l \mp\omega_l\ln (1\mp e^{-\alpha-\beta\epsilon_l})$$

当满足非简并条件时，$e^{-\alpha}\ll 1$，所以可以用泰勒展开：

$$\sum_l\mp\omega_l\ln (1\mp e^{-\alpha-\beta\epsilon_l}) \approx\sum_l\omega_l e^{-\alpha-\beta\epsilon_l}=e^{-\alpha}Z_1$$

可以看到，直接化为了玻尔兹曼统计的配分函数。

弱简并理想玻色气体和费米气体

在 玻尔兹曼统计 中，我们计算了理想气体的若干性质，现在我们考虑弱简并的理想玻色气体和费米气体。回顾此前的粒子数计算（用巨正则配分函数统一描述）：

$$N=-\frac{\partial }{\partial \alpha} \ln \Xi = -\frac{\partial e^{-\alpha}Z_1}{\partial \alpha}=e^{-\alpha}Z_1=\ln \Xi$$

这说明在非简并条件下，$\ln \Xi$对$\alpha$求导只改变符号。

在弱简并环境下就不能忽略：

$$\ln \Xi =\sum_l \mp\omega_l\ln (1\mp e^{-\alpha-\beta\epsilon_l})$$

的近似影响。代入 近独立粒子的最概然分布 中的自由粒子态密度公式，直接写出配分函数的计算公式：

$$\ln \Xi=\mp\int_0^\infty D(\epsilon)\ln (1\mp e^{-\alpha-\beta\epsilon})d\epsilon=\mp\frac{2\pi V}{h^3}(2m)^\frac32\int_0^\infty \epsilon^{1/2}\ln (1\mp e^{-\alpha-\beta\epsilon})d\epsilon$$

Note：对于玻色气体，如果温度太低，可能会出现$e^{-\alpha}\geq 2.612$的情况，此时玻色-爱因斯坦凝聚会发生，这时候上述式子就会因为$D(\epsilon=0)=0$而不成立。这一点在下一节会讨论。

考虑弱简并条件，我们可以对$\ln (1- e^{-\alpha-\beta\epsilon})$进行泰勒展开：

$$\ln (1- e^{-\alpha-\beta\epsilon})\approx - \sum_{n=1}^\infty \frac{e^{-n(\alpha+\beta\epsilon)}}{n}$$

代入上式，得到：

$$\begin{aligned}\ln \Xi &= \frac{2\pi V}{h^3}(2m)^\frac32\int_0^\infty \epsilon^{1/2}\sum_{n=1}^\infty \frac{e^{-n(\alpha+\beta\epsilon)}}{n}d\epsilon\\ &=\frac{2\pi V}{h^3}(2m)^\frac32\sum_{n=1}^\infty \frac{e^{-n\alpha}}{n}\int_0^\infty \epsilon^{1/2} e^{-n\beta\epsilon}d\epsilon\\ (x=n\beta\epsilon)&=\frac{V}{h^3}(\frac{2\pi m}{\beta})^\frac32\sum_{n=1}^\infty \frac{e^{-n\alpha}}{n^{5/2}}\\ &=Z_1g_{\frac52}(z)\\ \end{aligned}$$

这里引入理想气体的配分函数$Z_1$，逸度$z$和玻色函数$g_\nu(z)$：

$$g_\nu(z)=\frac{1}{\Gamma(\nu)}\int_0^\infty \frac{t^{\nu-1}}{e^t/z-1}dt=\sum_{n=1}^\infty \frac{z^n}{n^\nu}$$

现在，我们可以放心求解各个热力学量了：

$$\begin{cases} N&=-\dfrac{\partial z}{\partial \alpha} \dfrac{\partial \ln \Xi}{\partial z} =Z_1 g_{\frac32}(z)\\ U&=-\dfrac{\partial \ln \Xi}{\partial \beta} = \dfrac32 kT Z_1g_{\frac52}(z)\\ J&=-kT\ln \Xi = -kT Z_1 g_{\frac52}(z)\\ S&=k\left(\ln \Xi+\alpha N + \beta U\right)=k\left(\frac52Z_1 g_{\frac52}(z) + \alpha Z_1 g_{\frac32}(z)\right) \end{cases}$$

进行简单的变换，就是：

$$U=\frac32 NkT\frac{g_{\frac52}(z)}{g_{\frac32}(z)}$$ $$pV=-J=\frac32 U$$

同样的道理，可以定义费米函数：

$$f_\nu(z)=\frac{1}{\Gamma(\nu)}\int_0^\infty \frac{t^{\nu-1}}{e^t/z+1}dt=\sum_{n=1}^\infty \frac{(-1)^{n+1}z^n}{n^\nu}$$

导出费米气体的热力学量，这里不再赘述。显然，费米函数在同阶情况下小于玻色函数。

对比理想气体和理想玻色气体的热力学量：

量	理想气体	理想玻色气体	理想费米气体
配分函数	$Z_1=\dfrac{V}{h^3}\left(\dfrac{2\pi m}{\beta}\right)^\frac32$	$\ln\Xi=Z_1 g_{\frac52}(z)$	$\ln\Xi=Z_1 f_{\frac52}(z)$
$N$	$Z_1e^{-\alpha}=Z_1g_{\infty}(z)$	$Z_1 g_{\frac32}(z)$	$Z_1 f_{\frac32}(z)$
$U$	$\dfrac32 NkT$	$\dfrac32 NkT\dfrac{g_{\frac52}(z)}{g_{\frac32}(z)}$	$\dfrac32 NkT\dfrac{f_{\frac52}(z)}{f_{\frac32}(z)}$
$S$	$k\left(\dfrac52Z_1e^{-\alpha}+\alpha Z_1e^{-\alpha}\right)$	$k\left(\dfrac52Z_1 g_{\frac52}(z) + \alpha Z_1 g_{\frac32}(z)\right)$	$k\left(\dfrac52Z_1 f_{\frac52}(z) + \alpha Z_1 f_{\frac32}(z)\right)$

可以分析，理想玻色气体的内能较理想气体更小，而理想费米气体的内能较理想气体更大，这表明玻色子之间存在等效的吸引作用，而费米子之间存在等效的排斥作用。

import numpy as np
import matplotlib.pyplot as plt

def g_nu(z, nu, n_terms=1000):
    """计算玻色函数 g_nu(z) = sum_{n=1}^∞ z^n / n^nu (级数展开)"""
    n = np.arange(1, n_terms + 1)
    return np.sum(z**n / (n**nu))

def f_nu(z, nu, n_terms=1000):
    """计算费米函数 f_nu(z) = sum_{n=1}^∞ (-1)^{n+1} z^n / n^nu (级数展开)"""
    n = np.arange(1, n_terms + 1)
    return np.sum((-1)**(n + 1) * z**n / (n**nu))

# 定义 z 的范围 (0 < z < 1)
z = np.linspace(0.01, 0.99, 1000)

# 计算 g_{5/2}(z) / g_{3/2}(z)
g_ratio = np.array([g_nu(zi, 5/2) / g_nu(zi, 3/2) for zi in z])

# 计算 f_{5/2}(z) / f_{3/2}(z)
f_ratio = np.array([f_nu(zi, 5/2) / f_nu(zi, 3/2) for zi in z])

# 绘制图形
plt.figure(figsize=(10, 6))
plt.plot(z, g_ratio, label=r'$\frac{g_{5/2}(z)}{g_{3/2}(z)}$', color='blue')
plt.plot(z, f_ratio, label=r'$\frac{f_{5/2}(z)}{f_{3/2}(z)}$', color='red')
plt.axhline(y=1, color='black', linestyle='--', label='1')

plt.xlabel('z (fugacity)')
plt.ylabel('Ratio')
plt.title(r'Comparison of $\frac{g_{5/2}(z)}{g_{3/2}(z)}$ and $\frac{f_{5/2}(z)}{f_{3/2}(z)}$ with 1')
plt.legend()
plt.grid(True)
plt.show()

玻色-爱因斯坦凝聚(BEC)

由上述计算可知：

$$N=Z_1 g_{\frac32}(z)=\frac{V}{h^3}\left(\frac{2\pi m}{\beta}\right)^\frac32 g_{\frac32}(z)$$

显然，由于粒子数是守恒的，当温度降低的时候，$Z_1$会变小，而$g_{\frac32}(z)$会变大。遗憾的是，$g_{\frac32}(z)$是有上界$g_{\frac32}(z)\leq g_{\frac32}(1)=\zeta(\frac32)\approx 2.612$的，且当$z>1$的时候，玻色函数发散。这里的原因在上面也提到了，是因为计算配分函数的时候，$\epsilon=0$的态密度为0：

$$\ln \Xi=-\frac{2\pi V}{h^3}(2m)^\frac32\int_0^\infty \epsilon^{1/2}\ln (1- e^{-\alpha-\beta\epsilon})d\epsilon$$

这与玻色分布$a_l=\dfrac{\omega_l}{e^{\alpha+\beta\epsilon_l}-1}(e^{\alpha}\geq1,\mu\leq0)$相违背。

所以当温度降低到某个临界值时，$g_{\frac32}(z)$会达到上限，此时：

$$z=1,\alpha=0\Rightarrow \mu=0,T_c = \frac{2\pi\hbar^2}{mk}\left[\frac{n}{g_{\frac32}(1)}\right]^{2/3}$$

当温度继续降低，把配分函数的第一项提取出来：

$$\begin{aligned} \ln \Xi&=-\ln (1- e^{-\alpha})-\frac{2\pi V}{h^3}(2m)^\frac32\int_0^\infty \epsilon^{1/2}\ln (1- e^{-\alpha-\beta\epsilon})d\epsilon\\ &=-\ln (1- z)+Z_1 g_{\frac52}(z)\\ \end{aligned}$$

对其求导得到粒子数：

$$N=-\frac{\partial z}{\partial \alpha} \frac{\partial \ln \Xi}{\partial z} = \frac{z}{1-z}+Z_1 g_{\frac32}(z)=N_0+N_{\text{exc}}$$

这里的$N_0$是基态粒子数，$N_{\text{exc}}$是激发态粒子数。由于激发态粒子数目和温度成$N_{\text{exc}}\propto T^{3/2}$关系，这个式子可以转写为：

$$N=N_0+N_{\text{exc}}=N_0+N\left(\frac{T}{T_c}\right)^{3/2}$$

所以：

$$N_0(T)=N\left[1-\left(\frac{T}{T_c}\right)^{3/2}\right]$$

这时候激发态的能量为：

$$\begin{aligned} U&=-\frac{\partial \ln \Xi}{\partial \beta} \\ (z=1)&= -\ln (1- z)+\frac32 N_{\text{exc}}kT\frac{g_{\frac52}(z)}{g_{\frac32}(z)}\\ &= \frac32 NkT\left(\frac{T}{T_c}\right)^{3/2}\frac{g_{\frac52}(1)}{g_{\frac32}(1)}\\ &= 0.770NkT\left(\frac{T}{T_c}\right)^{3/2} \end{aligned}$$

综合大于临界温度的情况，内能表达式为：

$$U(T)=\begin{cases} \frac32 NkT\frac{g_{\frac52}(z)}{g_{\frac32}(z)} & T>T_c\\ 0.770NkT\left(\frac{T}{T_c}\right)^{3/2} & T<T_c \end{cases}$$

热容为（z也是T的函数，所以这里写不出解析的表达式）：

$$C_V=\frac{\partial U}{\partial T}=\begin{cases} \frac{\partial }{\partial T}\left(\frac{3}{2}NkT\frac{g_{\frac52}(z)}{g_{\frac32}(z)}\right) & T>T_c\\ 1.925Nk\left(\frac{T}{T_c}\right)^{1/2} & T<T_c \end{cases}$$

化学势为：

$$\mu = \begin{cases} 0 & T<T_c\\ kT\ln z & T\geq T_c \end{cases}$$

如果是二维的非相对论粒子，则积分变为：

$$N\propto \int \dfrac{1}{1-e^{-\beta\epsilon}}d\epsilon$$

这个积分是发散的，所以二维非相对论粒子不会出现玻色-爱因斯坦凝聚现象。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import zeta, gamma
from scipy.integrate import quad
from functools import lru_cache

# 物理常数
h = 6.626e-34       # 普朗克常数 (J·s)
m = 1.67e-27        # 氦原子质量 (kg)
k = 1.38e-23        # 玻尔兹曼常数 (J/K)
N = 1e24            # 总粒子数
V = 1e-3            # 体积 (m³)

# 预计算常数乘积
h2 = h**2
two_pi_mk = 2 * np.pi * m * k
prefactor_Tc = h2 / two_pi_mk
zeta_3_2 = zeta(3/2)

# 计算临界温度 T_c
n = N / V
T_c = prefactor_Tc * (n / zeta_3_2)**(2/3)

# 预计算用于玻色函数的常数
h3 = h**3
two_pi_mk_1_5 = (two_pi_mk)**1.5

# 玻色函数计算（积分形式）
@lru_cache(maxsize=1000)
def g_nu(z, nu):
    if z < 0.9:  # 对于较小的z，使用级数展开更高效
        return np.sum(z**np.arange(1, 1000) / np.arange(1, 1000)**nu)
    else:        # 对于接近1的z，使用积分形式
        def integrand(t):
            return t**(nu-1) / (np.exp(t)/z - 1)
        integral, _ = quad(integrand, 0, np.inf)
        return integral / gamma(nu)

# 缓存已计算的z值
z_cache = {}

# 二分法求解 z(T)，使用缓存
def solve_z(T):
    if T in z_cache:
        return z_cache[T]
    
    Z_1 = (V / h3) * two_pi_mk_1_5 * T**1.5
    left, right = 1e-16, 1 - 1e-16
    for _ in range(200):
        mid = (left + right)/2
        val = Z_1 * g_nu(mid, 1.5) - N
        if abs(val) < 1e-6: 
            z_cache[T] = mid
            return mid
        (left, right) = (left, mid) if val > 0 else (mid, right)
    
    z_cache[T] = (left + right)/2
    return z_cache[T]

# 内能计算，使用缓存
@lru_cache(maxsize=1000)
def U(T):
    if T > T_c:
        z = solve_z(T)
        g_1_5 = g_nu(z, 1.5)
        g_2_5 = g_nu(z, 2.5)
        return 1.5 * N * k * T * g_2_5 / g_1_5
    else:
        return 0.770 * N * k * T * (T/T_c)**1.5

# 数值差分计算热容，使用缓存
def C_V(T, delta_T=1e-4):
    if T > T_c + delta_T:
        return (U(T + delta_T) - U(T - delta_T)) / (2 * delta_T)
    elif T < T_c - delta_T:
        return (U(T + delta_T) - U(T - delta_T)) / (2 * delta_T)
    else:  # 临界区附近
        return (U(T + delta_T) - U(T)) / delta_T

# 化学势计算，使用缓存
@lru_cache(maxsize=1000)
def chemical_potential(T):
    if T <= T_c:
        return 0.0
    else:
        z = solve_z(T)
        return k * T * np.log(z)

# 温度范围
T_range = np.linspace(0.1 * T_c, 2 * T_c, 200)
U_values = [U(T) for T in T_range]
C_V_values = [C_V(T) for T in T_range]
mu_values = [chemical_potential(T) for T in T_range]

# 平滑处理热容曲线
def smooth_spikes(T_range, C_V_values, threshold=1.5):
    C_V_smooth = np.array(C_V_values.copy())
    n = len(C_V_smooth)
    
    derivatives = np.abs(np.diff(C_V_smooth))
    mean_derivative = np.mean(derivatives)
    
    spike_indices = []
    for i in range(1, n-1):
        if (np.abs(C_V_smooth[i+1] - C_V_smooth[i]) > threshold * mean_derivative and
            np.abs(C_V_smooth[i] - C_V_smooth[i-1]) > threshold * mean_derivative):
            spike_indices.append(i)
    
    for i in spike_indices:
        left_val = C_V_smooth[i-1] if i > 0 else C_V_smooth[i]
        right_val = C_V_smooth[i+1] if i < n-1 else C_V_smooth[i]
        C_V_smooth[i] = (left_val + right_val) / 2
    
    return C_V_smooth

C_V_smoothed = smooth_spikes(T_range, C_V_values)

# 绘图
plt.figure(figsize=(18, 6))

# 左子图：内能
plt.subplot(1, 3, 1)
plt.plot(T_range, U_values, 'b-', linewidth=2)
plt.axvline(x=T_c, color='r', linestyle='--', label=f'$T_c$ = {T_c:.2f} K')
plt.xlabel('Temperature $T$ (K)', fontsize=12)
plt.ylabel('Internal Energy $U(T)$ (J)', fontsize=12)
plt.title('Internal Energy of Ideal Bose Gas', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)

# 中子图：热容（平滑后）
plt.subplot(1, 3, 2)
plt.plot(T_range, C_V_values, 'g-', alpha=0.3, label='Original')
plt.plot(T_range, C_V_smoothed, 'b-', linewidth=2, label='Smoothed')
plt.axvline(x=T_c, color='r', linestyle='--', label=f'$T_c$ = {T_c:.2f} K')
plt.xlabel('Temperature $T$ (K)', fontsize=12)
plt.ylabel('Heat Capacity $C_V(T)$ (J/K)', fontsize=12)
plt.title('Smoothed Heat Capacity', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)

# 右子图：化学势
plt.subplot(1, 3, 3)
plt.plot(T_range, mu_values, 'b-', linewidth=2)
plt.axvline(x=T_c, color='r', linestyle='--', label=f'$T_c$ = {T_c:.2f} K')
plt.xlabel('Temperature $T$ (K)', fontsize=12)
plt.ylabel('Chemical Potential $\mu$ (J)', fontsize=12)
plt.title('Chemical Potential of Ideal Bose Gas', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

热辐射和光子气体

瑞利-金斯公式

在上面我们讨论的都是非相对论的能量关系，即$\epsilon=\frac{p^2}{2m}$。如果考虑相对论的能量关系（比如），我们就要改写 近独立粒子的最概然分布 中的粒子数密度公式：

$$\prod_i dn_i=\frac{4\pi Vp^2}{h^3} dp$$ $$\epsilon=cp\Rightarrow dp=d\epsilon/c$$ $$D(\epsilon)d\epsilon=\prod_i dn_i=\frac{4\pi V\epsilon^2}{c^3h^3} d\epsilon$$

一般我们改写为以频率为自变量：

$$D(\omega)d\omega = \frac{V}{2\pi^2c^3}\omega^2d\omega$$

考虑到光子具有自旋量子数为1，不存在$m_z=0$的光子，那么简并度$\omega_l=2$，所以光子的态密度为：

$$D(\omega)d\omega = \frac{V}{\pi^2c^3}\omega^2d\omega$$

瑞利-金斯根据能均分定理，认为每一个振动自由度具有能量$kT$，从而可以得到辐射场的能量谱：

$$\begin{cases}U(\omega,T)d\omega &= D(\omega)kTd\omega=\frac{V}{\pi^2c^3}\frac{\omega^2}{\beta}d\omega\\U(\lambda,T)d\lambda &= U(\frac{2\pi c}{\lambda},T)\frac{2\pi c}{\lambda^2}d\lambda=\frac{8\pi V}{\lambda^4}\frac{1}{\beta}d\lambda\end{cases}$$

瑞利近似在低频段拟合较好，但在高频段出现紫外灾难，即随着频率的增加，能量密度趋向无穷大，这显然不符合实际。归根结底还是将能均分定理用在所有自由度上是错误的。

维恩公式

维恩假设光子气体和理想气体一样，具有玻尔兹曼分布：

$$\begin{cases}U(\omega,T)d\omega &= D(\omega)\hbar\omega e^{-\beta\hbar\omega}d\omega = \frac{V}{\pi^2c^3}\hbar\omega^3e^{-\beta\hbar\omega}d\omega\\U(\lambda,T)d\lambda &= U(\frac{2\pi c}{\lambda},T)\frac{2\pi c}{\lambda^2}d\lambda=\frac{8\pi hc V}{\lambda^5}e^{-\frac{\beta hc}{\lambda}}d\lambda\end{cases}$$

维恩公式在高频段拟合较好，但在低频段出现了明显的偏差。

配分函数：

$$Z_1=\int_0^\infty D(\omega)d\omega e^{-\beta\hbar\omega}=\frac{V}{2\pi^2c^3}\int_0^\infty \omega^2e^{-\beta\hbar\omega}d\omega=\frac{V}{2\pi^2c^3}\frac{2}{\beta^3\hbar^3}$$

得到内能和压强：

$$\begin{cases}U(T) &= -\frac{\partial \ln Z_1}{\partial \beta} = \frac{3V}{2\pi^2c^3}\frac{2}{\beta^4\hbar^3}\\P(T) &= -\frac{\partial \ln Z_1}{\partial V} = \frac{}{2\pi^2c^3}\frac{2}{\beta^3\hbar^3}\end{cases}\Rightarrow U=3PV$$

对比理想气体的内能和压强：

$$U=\frac32 NkT=\frac32 PV$$

相对论的能量关系带来了系数的区别。

普朗克公式

使用光子气体描述热辐射场的观点是正确的，但是在确定光子气体的统计分布时，不能使用玻尔兹曼分布，更不能使用一致分布。光子气体的统计分布应该是玻色-爱因斯坦分布。由于光子气体是粒子数不守恒的，那么拉格朗日法只需要引入一个乘子：

$$\begin{cases}U(\omega,T)d\omega &= \frac{V}{\pi^2c^3}\frac{\hbar\omega^3d\omega}{e^{\beta\hbar\omega}-1}\\ U(\lambda,T)d\lambda &= \frac{8\pi hc V}{\lambda^5}\frac{1}{e^{\frac{\beta hc}{\lambda}}-1}d\lambda\end{cases}$$

$\alpha=-\frac{\mu}{kT}=0$意味着光子气体的化学势为0。可以看到，普朗克公式在低频段和高频段的近似分别就是瑞利-金斯公式和维恩公式。

总的内能为：

$$U(T) = \int_0^\infty U(\omega,T)d\omega = \frac{V}{\pi^2c^3}\int_0^\infty \frac{\hbar\omega^3d\omega}{e^{\beta\hbar\omega}-1}=\frac{\pi^2k^4}{15c^3\hbar^3}VT^4$$

内能密度的最大值有：

$$\frac{dU}{dT} = \frac{d}{dx}\frac{\omega^3d\omega}{e^{\beta\hbar\omega}-1}=0\Rightarrow \hbar\omega=2.822kT$$

这就是维恩位移定律。

也可以使用巨配分函数来计算，效果是一样的。光子气体的巨配分函数为：

$$\Xi = \prod_l (1-e^{-\beta\epsilon_l})^{-\omega_l} = \prod_l (1-e^{-\beta\hbar\omega})^{-2}$$ $$\ln \Xi = -2\int_0^\infty \ln(1-e^{-\beta\hbar\omega})d\omega=-\frac{V}{\pi^2c^3}\int_0^\infty \omega^2\ln(1-e^{-\beta\hbar\omega})d\omega=\frac{\pi^2V}{45c^3}\frac{1}{(\beta\hbar)^3}$$

所以内能：

$$U = -\frac{\partial J}{\partial \beta} = \frac{\pi^2}{15c^3}\frac{1}{\beta^4\hbar^3}$$

压强：

$$P = -\frac{\partial J}{\partial V} = \frac{1}{\beta}\frac{\partial \ln \Xi}{\partial V} = \frac{\pi^2}{45c^3}\frac{1}{\beta^4\hbar^3}$$

那么两者的关系：

$$P = \frac13 \frac UV$$

熵：

$$S = \frac{U+PV}{T} = \frac{4k\pi^2}{45c^3}\frac{V}{\beta^3\hbar^3}$$

金属中的自由电子气体

金属中的自由电子气体是强简并理想费米气体（这是因为费米能的量级远大于$kT$）。根据费米-狄拉克分布：

$$a_l = \frac{2}{e^{\beta(\epsilon-\mu)}+1}$$

这里考虑了自旋的两个简并。

电子的量子态数为：

$$D(\epsilon)d\epsilon = 4\pi V(\frac{2m}{\hbar^2})^{3/2}\epsilon^{1/2}d\epsilon$$

平均电子数为：

$$dn = 4\pi V(\frac{2m}{\hbar^2})^{3/2}\frac{\epsilon^{1/2}}{e^{\beta(\epsilon-\mu)}+1}d\epsilon$$

化学势通过粒子总数求得：

$$N = \int_0^\infty 4\pi V(\frac{2m}{\hbar^2})^{3/2}\frac{\epsilon^{1/2}}{e^{\beta(\epsilon-\mu)}+1}d\epsilon$$

对于0K的电子气：

$$\begin{cases}a_l=2,\epsilon<\mu(0)\\a_l=0,\epsilon>\mu(0)\end{cases}$$

所以：

$$N = \int_0^{\mu(0)} 4\pi V(\frac{2m}{\hbar^2})^{3/2}\epsilon^{1/2}d\epsilon$$

解得：

$$\mu(0) = \frac{\hbar^2}{2m}\left(\frac{3\pi^2N}{V}\right)^{2/3}$$

这就是费米能级，令$\mu(0)=\frac{p_F^2}{2m}$，得到费米动量：

$$p_F = \hbar\left(\frac{3\pi^2N}{V}\right)^{1/3}$$

对于任意温度的电子气，可以参考 金属电子论 。