第4章:固体の熱的性質

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# 物理定数
R = 8.314  # 気体定数 [J/(mol·K)]
k_B = 1.380649e-23  # Boltzmann定数 [J/K]
hbar = 1.054571817e-34  # Reduced Planck定数 [J·s]

class HeatCapacityModels:
    """固体の比熱を計算するモデルクラス"""

    def __init__(self, theta_E=None, theta_D=None):
        """
        Parameters:
        -----------
        theta_E : float
            Einstein温度 [K]
        theta_D : float
            Debye温度 [K]
        """
        self.theta_E = theta_E
        self.theta_D = theta_D

    def einstein_heat_capacity(self, T):
        """
        Einsteinモデルによる定容モル比熱を計算

        Parameters:
        -----------
        T : float or array
            温度 [K]

        Returns:
        --------
        C_V : float or array
            定容モル比熱 [J/(mol·K)]
        """
        if self.theta_E is None:
            raise ValueError("Einstein temperature not set")

        x = self.theta_E / T
        C_V = 3 * R * x**2 * np.exp(x) / (np.exp(x) - 1)**2
        return C_V

    def debye_integrand(self, x):
        """Debye積分の被積分関数"""
        if x < 1e-10:  # 0除算回避
            return 0
        return x**4 * np.exp(x) / (np.exp(x) - 1)**2

    def debye_heat_capacity(self, T):
        """
        Debyeモデルによる定容モル比熱を計算

        Parameters:
        -----------
        T : float or array
            温度 [K]

        Returns:
        --------
        C_V : float or array
            定容モル比熱 [J/(mol·K)]
        """
        if self.theta_D is None:
            raise ValueError("Debye temperature not set")

        # スカラー値の場合
        if np.isscalar(T):
            if T < 1e-10:
                return 0
            x_D = self.theta_D / T
            integral, _ = quad(self.debye_integrand, 0, x_D)
            C_V = 9 * R * (T / self.theta_D)**3 * integral
            return C_V

        # 配列の場合
        C_V = np.zeros_like(T)
        for i, temp in enumerate(T):
            if temp > 1e-10:
                x_D = self.theta_D / temp
                integral, _ = quad(self.debye_integrand, 0, x_D)
                C_V[i] = 9 * R * (temp / self.theta_D)**3 * integral
        return C_V

    def debye_T3_law(self, T):
        """
        Debye T³則（低温極限）

        Parameters:
        -----------
        T : float or array
            温度 [K]

        Returns:
        --------
        C_V : float or array
            定容モル比熱 [J/(mol·K)]
        """
        if self.theta_D is None:
            raise ValueError("Debye temperature not set")

        return (12 * np.pi**4 / 5) * R * (T / self.theta_D)**3


# 使用例: 銅の比熱
theta_D_Cu = 343  # 銅のDebye温度 [K]
theta_E_Cu = 240  # 銅の等価Einstein温度 [K]

model = HeatCapacityModels(theta_E=theta_E_Cu, theta_D=theta_D_Cu)

# 温度範囲を設定
T = np.linspace(1, 500, 200)

# 各モデルで計算
C_einstein = model.einstein_heat_capacity(T)
C_debye = model.debye_heat_capacity(T)
C_T3 = model.debye_T3_law(T)

# 可視化
plt.figure(figsize=(10, 6))
plt.plot(T, C_einstein, label='Einstein model', linewidth=2)
plt.plot(T, C_debye, label='Debye model', linewidth=2)
plt.plot(T, C_T3, '--', label='Debye T³ law', linewidth=2, alpha=0.7)
plt.axhline(y=3*R, color='gray', linestyle='--',
            label='Dulong-Petit (3R)', alpha=0.5)

plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('$C_V$ [J/(mol·K)]', fontsize=12)
plt.title('Heat Capacity Models for Copper', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.xlim(0, 500)
plt.ylim(0, 28)

plt.tight_layout()
plt.show()

# 低温領域のズームイン
T_low = np.linspace(1, 100, 100)

C_einstein_low = model.einstein_heat_capacity(T_low)
C_debye_low = model.debye_heat_capacity(T_low)
C_T3_low = model.debye_T3_law(T_low)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 線形スケール
ax1.plot(T_low, C_einstein_low, label='Einstein', linewidth=2)
ax1.plot(T_low, C_debye_low, label='Debye', linewidth=2)
ax1.plot(T_low, C_T3_low, '--', label='T³ law', linewidth=2)
ax1.set_xlabel('Temperature (K)', fontsize=12)
ax1.set_ylabel('$C_V$ [J/(mol·K)]', fontsize=12)
ax1.set_title('Low Temperature: Linear Scale', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 対数スケール
ax2.semilogy(T_low, C_einstein_low, label='Einstein', linewidth=2)
ax2.semilogy(T_low, C_debye_low, label='Debye', linewidth=2)
ax2.semilogy(T_low, C_T3_low, '--', label='T³ law', linewidth=2)
ax2.set_xlabel('Temperature (K)', fontsize=12)
ax2.set_ylabel('$C_V$ [J/(mol·K)]', fontsize=12)
ax2.set_title('Low Temperature: Log Scale', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.show()

# T³則の検証
T_verify = np.linspace(5, 50, 10)
C_debye_verify = model.debye_heat_capacity(T_verify)

# C_V / T³ をプロット（一定値になるはず）
plt.figure(figsize=(8, 5))
plt.plot(T_verify, C_debye_verify / T_verify**3, 'o-',
         markersize=8, linewidth=2)
plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('$C_V / T^3$ [J/(mol·K⁴)]', fontsize=12)
plt.title('Verification of T³ Law', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 代表的な物質のDebye温度
materials = {
    'Lead (Pb)': 105,
    'Copper (Cu)': 343,
    'Aluminum (Al)': 428,
    'Silicon (Si)': 645,
    'Diamond (C)': 2230
}

plt.figure(figsize=(12, 7))

for material, theta_D in materials.items():
    model_temp = HeatCapacityModels(theta_D=theta_D)
    T_range = np.linspace(1, 1000, 200)
    C_V = model_debye_heat_capacity(T_range)

    plt.plot(T_range/theta_D, C_V/(3*R), label=material, linewidth=2)

plt.axhline(y=1.0, color='gray', linestyle='--',
            label='Classical limit', alpha=0.5)
plt.xlabel('$T / \\Theta_D$ (reduced temperature)', fontsize=12)
plt.ylabel('$C_V / 3R$ (reduced heat capacity)', fontsize=12)
plt.title('Universal Debye Heat Capacity Curve',
          fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.xlim(0, 3)
plt.ylim(0, 1.1)

plt.tight_layout()
plt.show()