第4章: 相互作用系と相転移

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve k_B = 1.380649e-23 # J/K def mean_field_ising_equation(m, T, J, z, h, k_B): """平均場方程式""" if T == 0: return m - np.sign(J*z*m + h) beta = 1 / (k_B * T) return m - np.tanh(beta * (J * z * m + h)) def solve_magnetization(T, J, z, h, k_B, m_init=0.5): """自己無撞着方程式を解く""" sol = fsolve(mean_field_ising_equation, m_init, args=(T, J, z, h, k_B)) return sol[0] # パラメータ J = 1e-21 # J (相互作用強度) z = 4 # 2次元正方格子 T_c = J * z / k_B # 臨界温度 # 温度範囲 T_range = np.linspace(0.01, 2.5 * T_c, 200) # h = 0 での自発磁化 m_positive = [] m_negative = [] for T in T_range: # 正の解 m_pos = solve_magnetization(T, J, z, 0, k_B, m_init=0.9) m_positive.append(m_pos) # 負の解 m_neg = solve_magnetization(T, J, z, 0, k_B, m_init=-0.9) m_negative.append(m_neg) # 外部磁場ありの場合 h_values = [0, 0.1e-21, 0.3e-21, 0.5e-21] h_eV = [h / 1.602e-19 for h in h_values] colors = ['blue', 'green', 'orange', 'red'] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 自発磁化の温度依存性 ax1 = axes[0, 0] ax1.plot(T_range / T_c, m_positive, 'b-', linewidth=2, label='m > 0') ax1.plot(T_range / T_c, m_negative, 'r--', linewidth=2, label='m < 0') ax1.axvline(1, color='k', linestyle=':', linewidth=1, label='T_c') ax1.set_xlabel('T / T_c') ax1.set_ylabel('Magnetization m') ax1.set_title('自発磁化（h = 0）') ax1.legend() ax1.grid(True, alpha=0.3) # 外部磁場の効果 ax2 = axes[0, 1] for h, h_ev, color in zip(h_values, h_eV, colors): m_h = [solve_magnetization(T, J, z, h, k_B) for T in T_range] ax2.plot(T_range / T_c, m_h, color=color, linewidth=2, label=f'h = {h_ev:.2f} eV') ax2.axvline(1, color='k', linestyle=':', linewidth=1) ax2.set_xlabel('T / T_c') ax2.set_ylabel('Magnetization m') ax2.set_title('外部磁場による磁化') ax2.legend() ax2.grid(True, alpha=0.3) # 磁化率（h = 0での応答） ax3 = axes[1, 0] # χ = ∂m/∂h |_{h=0} epsilon = 1e-23 # 微小磁場 chi_values = [] for T in T_range: m0 = solve_magnetization(T, J, z, 0, k_B) m_eps = solve_magnetization(T, J, z, epsilon, k_B) chi = (m_eps - m0) / epsilon chi_values.append(chi) ax3.plot(T_range / T_c, np.array(chi_values) * 1e21, 'purple', linewidth=2) ax3.axvline(1, color='k', linestyle=':', linewidth=1, label='T_c') ax3.set_xlabel('T / T_c') ax3.set_ylabel('χ (10⁻²¹ J⁻¹)') ax3.set_title('磁化率（T_c で発散）') ax3.legend() ax3.grid(True, alpha=0.3) # 臨界領域の拡大 ax4 = axes[1, 1] T_critical = np.linspace(0.5 * T_c, 1.5 * T_c, 100) m_critical = [solve_magnetization(T, J, z, 0, k_B, m_init=0.9) for T in T_critical] ax4.plot(T_critical / T_c, m_critical, 'b-', linewidth=2) ax4.axvline(1, color='k', linestyle=':', linewidth=2, label='T_c') ax4.set_xlabel('T / T_c') ax4.set_ylabel('Magnetization m') ax4.set_title('臨界領域の磁化') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_ising_mean_field.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("=== Ising平均場理論 ===\n") print(f"相互作用強度 J = {J:.2e} J = {J/1.602e-19:.4f} eV") print(f"配位数 z = {z}") print(f"臨界温度 T_c = {T_c:.2f} K\n") # 異なる温度での磁化 test_temps = [0.5 * T_c, 0.9 * T_c, 1.0 * T_c, 1.5 * T_c] for T_test in test_temps: m_test = solve_magnetization(T_test, J, z, 0, k_B, m_init=0.9) print(f"T = {T_test:.2f} K (T/T_c = {T_test/T_c:.2f}): m = {m_test:.4f}") print("\n臨界指数:") print(" 磁化: m ~ (T_c - T)^β, β = 1/2 (平均場)") print(" 磁化率: χ ~ |T - T_c|^{-γ}, γ = 1 (平均場)")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve, curve_fit k_B = 1.380649e-23 def solve_m(T, J, z, h, k_B): """平均場方程式を解く""" if T < 1e-10: return np.sign(h) if h != 0 else 1.0 beta = 1 / (k_B * T) eq = lambda m: m - np.tanh(beta * (J * z * m + h)) return fsolve(eq, 0.5 if h >= 0 else -0.5)[0] J = 1e-21 z = 4 T_c = J * z / k_B # 臨界指数 β の計算 T_below_Tc = np.linspace(0.5 * T_c, 0.999 * T_c, 50) m_beta = [solve_m(T, J, z, 0, k_B) for T in T_below_Tc] reduced_temp_beta = (T_c - T_below_Tc) / T_c # べき乗フィッティング def power_law_beta(t, beta_exp): return t**beta_exp # 対数フィットで指数を求める log_t = np.log(reduced_temp_beta) log_m = np.log(np.abs(m_beta)) poly_beta = np.polyfit(log_t, log_m, 1) beta_fitted = poly_beta[0] # 臨界指数 γ の計算（磁化率） T_chi = np.linspace(1.01 * T_c, 2 * T_c, 50) epsilon_h = 1e-23 chi_gamma = [] for T in T_chi: m0 = solve_m(T, J, z, 0, k_B) m_h = solve_m(T, J, z, epsilon_h, k_B) chi = (m_h - m0) / epsilon_h chi_gamma.append(chi) reduced_temp_gamma = (T_chi - T_c) / T_c log_t_gamma = np.log(reduced_temp_gamma) log_chi = np.log(chi_gamma) poly_gamma = np.polyfit(log_t_gamma, log_chi, 1) gamma_fitted = -poly_gamma[0] # 臨界指数 δ の計算（T = T_c） h_range = np.logspace(-24, -20, 30) m_delta = [solve_m(T_c, J, z, h, k_B) for h in h_range] log_h = np.log(h_range) log_m_delta = np.log(np.abs(m_delta)) poly_delta = np.polyfit(log_h, log_m_delta, 1) delta_fitted = 1 / poly_delta[0] # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 臨界指数 β ax1 = axes[0, 0] ax1.loglog(reduced_temp_beta, m_beta, 'bo', markersize=6, label='Data') ax1.loglog(reduced_temp_beta, reduced_temp_beta**0.5, 'r--', linewidth=2, label=f'β = 0.5 (theory)') ax1.loglog(reduced_temp_beta, reduced_temp_beta**beta_fitted, 'g-', linewidth=2, label=f'β = {beta_fitted:.3f} (fit)') ax1.set_xlabel('(T_c - T) / T_c') ax1.set_ylabel('Magnetization m') ax1.set_title('臨界指数 β') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # 臨界指数 γ ax2 = axes[0, 1] ax2.loglog(reduced_temp_gamma, np.array(chi_gamma) * 1e21, 'go', markersize=6, label='Data') ax2.loglog(reduced_temp_gamma, 1e21 * (reduced_temp_gamma)**(-1), 'r--', linewidth=2, label='γ = 1 (theory)') ax2.loglog(reduced_temp_gamma, 1e21 * (reduced_temp_gamma)**(-gamma_fitted), 'b-', linewidth=2, label=f'γ = {gamma_fitted:.3f} (fit)') ax2.set_xlabel('(T - T_c) / T_c') ax2.set_ylabel('χ (10⁻²¹ J⁻¹)') ax2.set_title('臨界指数 γ') ax2.legend() ax2.grid(True, alpha=0.3, which='both') # 臨界指数 δ ax3 = axes[1, 0] ax3.loglog(h_range / 1.602e-19, np.abs(m_delta), 'mo', markersize=6, label='Data') ax3.loglog(h_range / 1.602e-19, (h_range/1e-21)**(1/3), 'r--', linewidth=2, label='δ = 3 (theory)') ax3.loglog(h_range / 1.602e-19, (h_range/1e-21)**(1/delta_fitted), 'c-', linewidth=2, label=f'δ = {delta_fitted:.3f} (fit)') ax3.set_xlabel('h (eV)') ax3.set_ylabel('m (at T = T_c)') ax3.set_title('臨界指数 δ') ax3.legend() ax3.grid(True, alpha=0.3, which='both') # 比熱（数値微分） ax4 = axes[1, 1] T_heat = np.linspace(0.5 * T_c, 1.5 * T_c, 100) dT = T_heat[1] - T_heat[0] free_energy = [] for T in T_heat: m_eq = solve_m(T, J, z, 0, k_B) # 自由エネルギー（平均場） if T > 0: beta = 1 / (k_B * T) f = -J * z * m_eq**2 / 2 - k_B * T * np.log(2 * np.cosh(beta * J * z * m_eq)) else: f = -J * z free_energy.append(f) entropy = -np.gradient(free_energy, dT) heat_capacity = T_heat * np.gradient(entropy, dT) ax4.plot(T_heat / T_c, -heat_capacity / k_B, 'r-', linewidth=2) ax4.axvline(1, color='k', linestyle=':', linewidth=2, label='T_c') ax4.set_xlabel('T / T_c') ax4.set_ylabel('C / k_B') ax4.set_title('比熱（T_c で不連続）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_critical_exponents.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("=== 臨界指数（平均場理論）===\n") print(f"理論値（平均場）:") print(f" β = 1/2 （磁化）") print(f" γ = 1 （磁化率）") print(f" δ = 3 （臨界等温線）") print(f" α = 0 （比熱、不連続）\n") print(f"数値フィット:") print(f" β = {beta_fitted:.4f}") print(f" γ = {gamma_fitted:.4f}") print(f" δ = {delta_fitted:.4f}\n") print("実験値（3次元Ising）:") print(" β ≈ 0.325") print(" γ ≈ 1.24") print(" δ ≈ 4.8") print(" α ≈ 0.11")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize_scalar def landau_free_energy(m, T, T_c, a0, b, h): """Landau自由エネルギー""" a = a0 * (T - T_c) return a * m**2 + b * m**4 - h * m def equilibrium_magnetization(T, T_c, a0, b, h): """平衡磁化（自由エネルギー最小化）""" # 複数の初期値から探索 m_candidates = [] for m_init in [-1, 0, 1]: result = minimize_scalar(lambda m: landau_free_energy(m, T, T_c, a0, b, h), bounds=(-2, 2), method='bounded') m_candidates.append(result.x) # 最小の自由エネルギーを与える解を選択 F_values = [landau_free_energy(m, T, T_c, a0, b, h) for m in m_candidates] return m_candidates[np.argmin(F_values)] # パラメータ T_c = 500 # K a0 = 1e-4 b = 1e-6 h_values = [0, 0.01, 0.05, 0.1] T_range = np.linspace(300, 700, 100) colors = ['blue', 'green', 'orange', 'red'] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Landauポテンシャルの形状（異なる温度） ax1 = axes[0, 0] m_plot = np.linspace(-1.5, 1.5, 200) temperatures_plot = [0.8*T_c, 0.95*T_c, T_c, 1.2*T_c] for T_val, color in zip(temperatures_plot, colors): F_plot = [landau_free_energy(m, T_val, T_c, a0, b, 0) for m in m_plot] ax1.plot(m_plot, F_plot, color=color, linewidth=2, label=f'T/T_c = {T_val/T_c:.2f}') ax1.set_xlabel('Order parameter m') ax1.set_ylabel('Free energy F(m)') ax1.set_title('Landauポテンシャル（h = 0）') ax1.legend() ax1.grid(True, alpha=0.3) # 平衡磁化の温度依存性 ax2 = axes[0, 1] for h, color in zip(h_values, colors): m_eq = [equilibrium_magnetization(T, T_c, a0, b, h) for T in T_range] ax2.plot(T_range / T_c, m_eq, color=color, linewidth=2, label=f'h = {h:.2f}') ax2.axvline(1, color='k', linestyle=':', linewidth=1, label='T_c') ax2.set_xlabel('T / T_c') ax2.set_ylabel('Magnetization m') ax2.set_title('平衡磁化') ax2.legend() ax2.grid(True, alpha=0.3) # 1次相転移（bを負にする場合） ax3 = axes[1, 0] # F = a m² + b m⁴ + c m⁶ (b < 0の場合、1次転移) b_first_order = -2e-6 c = 1e-7 def landau_first_order(m, T, T_c, a0, b, c, h): a = a0 * (T - T_c) return a * m**2 + b * m**4 + c * m**6 - h * m m_1st = np.linspace(-1.5, 1.5, 200) T_1st = [0.95*T_c, 0.98*T_c, T_c, 1.02*T_c] for T_val, color in zip(T_1st, colors): F_1st = [landau_first_order(m, T_val, T_c, a0, b_first_order, c, 0) for m in m_1st] ax3.plot(m_1st, F_1st, color=color, linewidth=2, label=f'T/T_c = {T_val/T_c:.2f}') ax3.set_xlabel('Order parameter m') ax3.set_ylabel('Free energy F(m)') ax3.set_title('1次相転移（二重井戸ポテンシャル）') ax3.legend() ax3.grid(True, alpha=0.3) # 磁化率 ax4 = axes[1, 1] epsilon_h = 0.001 chi_landau = [] for T in T_range: m0 = equilibrium_magnetization(T, T_c, a0, b, 0) m_h = equilibrium_magnetization(T, T_c, a0, b, epsilon_h) chi = (m_h - m0) / epsilon_h chi_landau.append(chi) ax4.plot(T_range / T_c, chi_landau, 'purple', linewidth=2) ax4.axvline(1, color='k', linestyle=':', linewidth=1, label='T_c') ax4.set_xlabel('T / T_c') ax4.set_ylabel('χ') ax4.set_title('磁化率（T_c で発散）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_landau_theory.png', dpi=300, bbox_inches='tight') plt.show() print("=== Landau理論 ===\n") print(f"パラメータ:") print(f" T_c = {T_c} K") print(f" a₀ = {a0}") print(f" b = {b}\n") print("2次相転移（b > 0）:") print(f" T > T_c: m = 0（単一井戸）") print(f" T < T_c: m = ±√(-a/2b)（二重井戸）") print(f" 臨界指数 β = 1/2\n") print("1次相転移（b < 0, c > 0）:") print(f" 二重井戸ポテンシャル") print(f" 転移で磁化が不連続にジャンプ") print(f" 潜熱が存在")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve k_B = 1.380649e-23 def vdw_pressure(V, T, a, b, k_B): """van der Waals圧力""" return k_B * T / (V - b) - a / V**2 def vdw_equation(V, P, T, a, b, k_B): """van der Waals方程式（体積を求める）""" return P - vdw_pressure(V, T, a, b, k_B) # van der Waalsパラメータ（Arガス） a = 1.355e-49 # J·m³ b = 3.201e-29 # m³ T_c = 8 * a / (27 * k_B * b) P_c = a / (27 * b**2) V_c = 3 * b print(f"臨界定数:") print(f" T_c = {T_c:.2f} K") print(f" P_c = {P_c/1e6:.2f} MPa") print(f" V_c = {V_c*1e30:.2f} ų\n") # 等温線（換算変数） v_range = np.linspace(0.4, 5, 300) temperatures_reduced = [0.85, 0.95, 1.0, 1.1, 1.3] colors = plt.cm.viridis(np.linspace(0, 1, len(temperatures_reduced))) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # P-V図（換算変数） ax1 = axes[0, 0] for t_r, color in zip(temperatures_reduced, colors): p_values = [(3 / v**2) * (8*t_r / (3*v - 1) - 1) for v in v_range if 3*v > 1] v_plot = [v for v in v_range if 3*v > 1] ax1.plot(v_plot, p_values, color=color, linewidth=2, label=f't = {t_r:.2f}') ax1.axhline(1, color='k', linestyle=':', linewidth=1) ax1.axvline(1, color='k', linestyle=':', linewidth=1) ax1.plot(1, 1, 'ro', markersize=10, label='Critical point') ax1.set_xlabel('v = V/V_c') ax1.set_ylabel('p = P/P_c') ax1.set_title('van der Waals等温線（換算変数）') ax1.set_xlim([0.4, 3]) ax1.set_ylim([0, 2]) ax1.legend() ax1.grid(True, alpha=0.3) # Maxwellの等面積則 ax2 = axes[0, 1] t_coexist = 0.90 v_fine = np.linspace(0.4, 3, 500) p_fine = [(3 / v**2) * (8*t_coexist / (3*v - 1) - 1) for v in v_fine if 3*v > 1] v_plot_fine = [v for v in v_fine if 3*v > 1] ax2.plot(v_plot_fine, p_fine, 'b-', linewidth=2, label=f't = {t_coexist}') # 平衡圧力の推定（簡易版） p_equilibrium = 0.7 # 目視で調整 ax2.axhline(p_equilibrium, color='r', linestyle='--', linewidth=2, label='Maxwell construction') ax2.set_xlabel('v = V/V_c') ax2.set_ylabel('p = P/P_c') ax2.set_title('Maxwell等面積則') ax2.set_xlim([0.4, 3]) ax2.set_ylim([0.4, 1.2]) ax2.legend() ax2.grid(True, alpha=0.3) # 圧縮率 ax3 = axes[1, 0] for t_r, color in zip(temperatures_reduced, colors): # κ = -1/V (∂V/∂P)_T v_comp = np.linspace(0.5, 3, 100) dv = v_comp[1] - v_comp[0] p_comp = [(3 / v**2) * (8*t_r / (3*v - 1) - 1) for v in v_comp if 3*v > 1] v_comp_valid = [v for v in v_comp if 3*v > 1] if len(p_comp) > 2: dP_dV = np.gradient(p_comp, dv) kappa = -1 / (np.array(v_comp_valid) * np.array(dP_dV)) ax3.plot(v_comp_valid, kappa, color=color, linewidth=2, label=f't = {t_r:.2f}') ax3.set_xlabel('v = V/V_c') ax3.set_ylabel('κ (reduced)') ax3.set_title('等温圧縮率') ax3.set_ylim([0, 10]) ax3.legend() ax3.grid(True, alpha=0.3) # 相図（T-P平面） ax4 = axes[1, 1] # 蒸気圧曲線（簡易推定） T_sat = np.linspace(0.6 * T_c, T_c, 50) # Clausius-Clapeyron近似 P_sat = P_c * np.exp(-1.5 * (1 - T_sat / T_c)) ax4.plot(T_sat / T_c, P_sat / P_c, 'b-', linewidth=3, label='Coexistence curve') ax4.plot(1, 1, 'ro', markersize=12, label='Critical point') ax4.fill_between(T_sat / T_c, 0, P_sat / P_c, alpha=0.3, color='blue', label='Liquid') ax4.fill_between(T_sat / T_c, P_sat / P_c, 2, alpha=0.3, color='red', label='Gas') ax4.set_xlabel('T / T_c') ax4.set_ylabel('P / P_c') ax4.set_title('相図（液-気共存曲線）') ax4.set_xlim([0.6, 1.2]) ax4.set_ylim([0, 1.5]) ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_vdw_phase_transition.png', dpi=300, bbox_inches='tight') plt.show() print("van der Waals気体の特徴:") print(" - 臨界点以下で液-気共存") print(" - Maxwell等面積則で平衡圧力決定") print(" - 臨界点で等温圧縮率が発散") print(" - 格子気体モデルと数学的に等価")